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PREFACE 


The  following  report  represents  the  third  review  of 
research  conducted  under  the  auspices  of  the  Joint  Services 
Electronics  Program  at  the  Institute  for  Electronics  Science 
at  Texas  Tech  University.  Specific  topics  covered  include, 
fault  analysis,  large-scale  systems,  stochastic  control  and 
estimation,  nonlinear  control,  multidimensional  system  theory, 

i 

optical  noise,  and  pattern  recognition. 
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ABSTRACT 

The  fault  diagnosis  problem  for  a linear  system  whose  transfer  function 
matrix  Is  measursd  at  a dlserets  sat  of  frequencies  Is  formalized.  A measure 
of  solvability  for  the  resultant  equations  and  a measure  of  testability  for 
the  unit  under  test  la  developed.  These,  In  turn,  are  used  as  the  basis  of 
algorithms  for  choosing  teat  points  and  test  frequencies. 

INTRODUCTION 

Conceptually,  the  fault  analysis  problem  for  an  analog  circuit  or  system 
amounts  to  the  measurement  of  a set  of  externally  accessible  parameters  of 
the  system  from  which  one  desires  to  determine  the  Internal  system  parameters 
or  equivalently  locate  Che  failed  components  as  illustrated  In  Figure  1. 


pi  r» 

1 2 3 


n-1  . • 


Figure  1.  Conceptual  Model  of  Fault  Diagnosis  Problem. 


Here,  the  measu 


ta,  m. 


represent  data  taken  at  distinct  test  points 


or  alternatively,  data  taken  at  a fixed  test  point  under  different  stimuli. 
Similarly,  the  ^ represent  parameters  characterizing  the  various  internal 
system  components.  Here,  a single  parameter  may  characterize  an  entire  component, 
a*T  * resistance,  capacitance  or  inductance.  Alternatively,  a component  may  be 


1 


2 


represented  by  several  parameters:  the  h-paramaters  of  a transistor,  the 
poles  and  gain  of  an  op— amp , etc*  In  general,  one  models  a system  component 
by  the  minimum  number  of  parameters  which  will  allow  the  failure  to  be  Iso- 
lated up  to  a shop  replaceable  assembly  (SRA)  with  all  "allowed"  system  fail- 
ures manifesting  themselves  In  the  form  of  some  parameter  change. 

To  solve  the  fault  diagnosis  problem,  one  then  measures  a * col(m^)  and 
solves  a nonlinear  algebraic  equation 

1.  m • F(r) 

for  r ■ col(r^)  to  diagnose  the  fault.  The  parameters  in  the  resultant  r 
vector  which  are  out  of  tolerance  then  indicate  the  faulty  component. ^ 

The  purpose  of  the  present  paper  is  to  give  an  explicit  formulation  of 
Che  fault  diagnosis  equations  which  arise  in  the  malntanence  of  linear  systems. 
Here,  one  measures  the  system  frequency  response  as  observed  from  a specified 
sec  of  externally  accessible  test  points  at  a discrete  sec  of  frequencies 
and  it  is  desired  to  solve  for  a vector  of  internal  system  parameters,  r, 
which  completely  characterize  the  frequency  response  matrices  of  the  in- 
dividual system  components;  Z^(s,r),  i • 1,  2,  ...,  q. 

In  the  following  section  the  explicit  form  for  the  fault  diagnosis 
equations  is  derived  for  a given  set  of  test  frequencies.  A measure  of  solva- 
bilitr  of  these  equations  is  then  developed  in  sectlon3  and  empolyed  in 
section  4 in  an  algorithm  for  optimally  selecting  test  frequencies.  The  measure 

of  solvability  for  the  fsnlt  analysis  equations , given  an  optimal  choice  of 

12  5 

cast  frequencies,  is  than  taken  as  a measure  of  testability  * * for  the  unit 

under  test  CUUT)  and  is  used  as  the  basis  of  an  algorithm  for  the  optical  choice  of 
3 4 5 

test  points.  ' ' Finally,  a number  of  examples  are  presented  in  seelon  5. 
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EXPLICIT  FORM  OF  THE  FAULT  DIAGNOSIS  EQUATIONS 


In  Che  ease  of  a linear  time- invariant  circuit  or  system,  the  fault 
diagnosis  equations  nay  be  expressed  In  analytical  form.^  Since  the  fault 
diagnosis  equations  deal  with  the  relationship  between  the  externally  measu 


able  ays  ten  parameters,  a,  and  the  Internal  component  parameters,  r,  we 

adopt  a component  connection  model  as  the  starting  point  for  the  derivation 

7 8 

of  Che  fault  diagnosis  equations.  * This  is  one  of  several  conmonly  em- 
ployed large  scale  system  models  In  which  the  conponents  and  connections  in  a 
circuit  or  system  are  modeled  by  distinct  equations,  thereby  permitting  one  to 
exp 11 cl tally  deal  with  Che  relationship  between  the  Individual  component 
parameters  and  the  composite  systen  parameters. 

Since  the  present  study  Is  restricted  to  linear  time- invariant  systems, 
we  assume  that  each  component  is  characterized  by  a transfer  function  matrix 
which  if  dependent  on  Che  potentially  variable  component  parameters,  Z^(s,r). 
For  the  classical  RLC  components  Z^(s,r)  may  take  the  form  R,  La,  or  1/sC 
for  the  case  of  a resistor,  inductor,  or  capacitor,  respectively.  More 
generally,  one  may  model  an  op-amp  by  the  transfer  function  k/(s-p1> (s-p2> 
where  the  parameter  vector,  r,  now  represents  the  three  potentially  variable 

mT 

component  parameters;  k,  p^,  p^;  or  a delay  by  ka  , etc.  Although  the  symbol 
Z Is  used,  the  components  are  not  assumed  to  be  represented  by  Impedance 


matrices.  Indeed,  hybrid  models  are  used  in  most  of  our  examples.  For  the 
purpose  of  analysis.  It  is  assumed  that  all  faults  manifest  themselves  la 
the  form  of  changes,  possibly  eataatophlc.  In  the  parameter  vector,  r,  with 


the  frequency  characteristics  of  the  components  unchanged.  Although  not 
universal,  this  fault  hypothesis  covers  the  most  commonly  encountered  situ- 
ations and  subsumes  the  common  industrial  practice  of  assuming  that  all 

failures  in  analog  circuits  and  systems  take  the  form  of  open  and  short 

a 

circuited  components. 
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Our  system  components  are  thus  characterized  by  a set  of  simultaneous 
equations 

2.  - Z^(s,r)a^  1 - 1,  2,  ....  q 

where  a^  and  bt  denote  the  conponent  Input  and  output  vectors , respectively. 

For  notation al  brevity*  these  component  equations  may  be  coofclned  into  a single 
block  diagonal  matrix  equation 

3.  b - Z(s ,r) a 

where  b ■ col(b^) , a » col(a^)  and  Z(s,r)  ■ dlag  (Z^(s,r)). 

Al chough  there  are  many  ways  to  represent  Che  connection  in  a circuit  or 
system;  say  a block  diagram,  linear  graph  or  signal  flow  graph,  any  such  repre- 
sentation is  sl^>ly  a graphical  means  for  displaying  a set  of  connection  e- 
a nations:  Kirchodf  lars,  adder  equations,  etc.  As  such,  for  our  conponent 
connection  model  we  adopt  a purely  algebraic  connection  model  in  which  the 
connection  equations  are  displayed  exp  lid  tally  without  the  intermediary  of 
some  kind  of  graphical  connection  diagram.  This  takes  the  form 

*.  * ■ Hi”  + Hz” 

■>  • Hi”  * H20 

where  u and  y represent  the  vectors  of  accessible  Inputs  and  outputs  which  are 

available  to  the  test  system.  In  slop  la  systems,  the  connection  matrices,  , 

an  usually  obtainable  by  inspection,  whenas.  In  non  coup  lex  systems, 

computer  codas  have  been  developed  for  their  derivation. ^ Moreover,  they  an 

8 

assund  to  exist  In  all  but  Che  most  pathaloglcal  systems. 

It  is  the  pair  of  simultaneous  matrix  equations  3 and  4 which  an  termed 
the  component  connection  model.  By  combining  equations  3 and  4 to  eliminate 
che  component  input  and  output  variables,  a and  b,  one  may  derive  * an  expnssion 


for  the  transfer  fraction  matrix  observable  by  the  test  system  between  the  test 


t 
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input  and  output  vectors,  u and  y,  obtaining 
t 5.  S(s,r)  - 42  + LjjU  - Z(s , r) L^) -1 2(s , r) 

where 

6.  y » S(a,r)u 

Por  a linear  tin*- invariant  system  the  transfer  function  S(s,r)  is  a 
couplets  description  of  the  mensurable  data  about  the  unit  under  test 
available  to  the  test  STSteo.  Moreover,  being  rational  it  is  completely  de- 
termined by  its  value  at  a finite  number  of  frequencies.  As  such,  without  loss 
of  generality,  we  may  take  our  measured  data  to  be  of  the  form 

7.  col[S(s^,r),  S(s2,r),  ...,  S(s^,r)] 

The  fault  diagnosis  equations  Chen  take  the  form 

Lj2  * l21(l-Z(.l,r)L11)"12(a1,c)L12 
42  + 41U-Z(>2.')t.11)'12(,2.r)L12 


SCsj^.r) 

S(s2,r) 


hi  + 41<l-2<«k.r)tu)"lZ(sktr>l1J 

Since  S(s, r)  is,  in  general,  a matrix,  the  fault  diagnosis  equations  as 
derived  above  take  the  form  of  a matrix  (col[S(s^,r) ] )valued  function  of  a vector 
valued  variable,  r.  Computationally , however,  we  prefer  to  work  with  a vector 
valued  function  of  a vector  valued  variable  and  hence,  we  transform  S(s,r)  into 
a co luma  vector  via 


SCs^.r) 


» 


r 
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9.  vec[S(s,r)]  - Col  [S^Cs.r)] 

where  S^(a,r)  denotes  the  ith  column  of  the  matrix,  S(s,r).  With  the  aid  of 

7 1 ^ 

the  identity  vec[XTZ]  -[Z%X]vec  [Y]  equation  8.  then  transforms  into  * 

- vecfLjj]  + [L^  ® L21(l-Z(s1,r)Lu)“1]vec[Z(s1,r)] 
• vecfL^]  + CL12  ® I*2^(l"Z(s2 »r)L^^)  3vecCZ(s2 » t)  ] 

e 

e 

- vecCL^  + 8 L2i(1“Z(V)Lll)~1^  v*cCZ(ak,r>^ 

which  is  the  form  of  the  fault  diagnosis  equations  with  which  we  desire  to  work. 

SOLVABILITY  OF  THE  FAULT  DIAGNOSIS  EQUATIONS 

For  the  fault  diagnosis  equations  derived  above  to  be  a viable  tool  of 
circuit  and  system  diagnosis  two  fundamental  questions  remain  to  be  answered: 
"What  test  frequencies  should  be  employed  to  optimize  the  solvability  of  the 
equations? "and  "How  solvable  are  the  equations  given  an  optimal  choice  of  test 
frequencies?"  Both  of  these  questions,  in  turn,  hinge  on  the  development  of  some 
type  of  measure  of  solvability  for  the  fault  diagnosis  equations. 

For  a sat  of  Linear  equations 

U.  m - Fr 

where  r is  an  n-vector,  m is  a p-vector  and  F is  a p by  a matrix  one  may 
characterize  the  solvability  of  the  equations  in  terms  of  the  number  of  arbitrary 
parameters  in  its  solution  (if  a solution  exists).  As  such,  5 • n-rank(F)  is  a 
natural  measure  of  the  solvability  for  equation  11.  Here,  5*0  implies  that  the 
equation  ham  a unique  solution,  <S  ■ 1 lap  lies  that  the  solution  is  determined  up 

i.  Ml  I 


vec[S(s1#r)] 

vec[S(s2,r)] 

e 

e 

e 

vec[S(sk,r)] 


- F(r) 
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to  one  arbitrary  par  am  car  ad  so  on,  with  Increasing  values  of  3 representing  de- 
creasing degrees  of  solvability. 

Uhfortwately,  the  fault  diagnosis  equations  are  nonlinear  even  for  linear 
system  ad  hence  we  oust  resort  to  the  Implicit  function  theorem  to  obtain  a 

jj 

manure  of  solvability  analogous  to  the  above.  Indeed,  If  r^  Is  a solution  to 
the  fault  diagnosis  equations,  then  r ^ la  determined  up  to  a 

12.  «(rf)  - n - rsnkjj—  (rf)]J 

dimensional  manifold  (of  arbitrary  parameters)  in  a neighborhood  of  rf.  Here 
dF/dr  is  the  Jacobla  matrix  of  partial  derivatives  of  F with  respect  to  r.  Vlth 
Che  aid  of  the  matrix  Identity  d(M  ^)/dr  ■ -M  ^[dM/dr]M  \ dF/dr  can  be  computed 


8 


In  turn,  allows  us  Co  tr ans fora  ths  local  measure  of  solvability  of  equation  12. 

into  a global  measure  of  solvability.  For  Chls  purpose  we  adopt  the  algebraic 

geometric  definition  for  Che  term  "almost  constant."  l.e.  we  say  that  a function 

of  r^  is  alnoat  constant  if  It  is  constant  except  possibly  for  those  values  of 

lying  in  an  algebraic  variety  (the  solution  specs  of  a finite  set  of  non-zero 

simultaneous  polynomial  equations  In  n variables).  More  generally*  we  say  that  a 

property  holds  "almost  everywhere"  (a.e. ) or  for  almost  all  r ^ In  n-spaca  If  It  is 

true  for  all  values  of  r^  except  possibly  those  lying  in  an  algebraic  variety. 

Since  the  Lebesque  measure  of  an  algebraic  variety  is  zero*  this  definition  for  the 

concept  "almost  everywhere"  is  consistent  with  the  more  common  measure  theoretic 

14 

definition  and  is  more  natural  in  the  context  of  our  application. 

Theorem  1:  Let  Z^s.r);  1 - 1,  2,  . ..,  q;  be  rational  in  r.  Then  6 (rf)  is  almost 
constant. 

Note,  che  assumption  that  Z^(s,r)  is  rational  in  r is  quite  minor  being  satis- 
fied by  all  of  the  examples  given  in  section  IX  except  for  the  delay  (which  can 
be  approximated  by  a function  which  Is  rational  in  r).  In  practice*  che 
component  transfer  function  matrices  will  also  be  rational  in  s though  chls  is  not 
required  for  che  present  theorem  since  F and  dF/dr  are  formulated  in  terms  of 

specific  test  frequencies*  s^*  » •••*  s^.  Given  our  assumption  on  the  Z^(s*r), 

dF 

together  with  equation  13.,  it  then  follows  Chat  -^(r^)  is  also  rational  in  r^. 
Proof  of  Theorem  1:  We  begin  by  showing  chat  an  arbitrary  polynomial  matrix  in 
r,  P(r) , has  almost  constant  rank.  Since  rank  P(r)  is  restricted  to  che  finite 
set  of  integers  (0,  1,  2,  ...,  j;  where  J is  the  minimum  of  the  ntaaber  of  ran 
and  eolums  in  P(r) ) * chare  exists  an  r which  maximizes  che  rank  of  P(r) 

B 

14.  ' rank[P(r  )]  > rankCP(r)] 

B 

Now,  the  rank  of  a matrix  is  the  dimension  of  its  largest  non-singular  square 
sub-matrix.  As  such,  P(r)  admits  a square  sub-matrix,  M(r),  whose  dimension  is 


9 


equal  co  che  rank  P(r J and  for  which 


' dee  M(r  ) »»  0. 


Now,  dec[M(r)]  la  a polynomial  in  r which  is  not  i den  el  c ally  zero  (from  equadon 

15.)  and  hence*  it  is  non- zero  a.e.  As  such. 


16.  rank[P(ra|)]  _>  rank[P(r)]  >,  rank[M(r)]  - rank[P(ra)]  a.e. 

showing  chat  rank[P(r)]  - rankfPCr^)]  almost  everywhere.  As  such,  rank[P(r)]is 
almost  constant.  y 

JP 

Sow,  to  verify  that  rank  [-^( r^)  ] is  constant  we  decompose  this  matrix  as 

17.  dP,_  „ P(rf) 

" d(rf) 

where  P(r^)  is  a polynomial  matrix  and  d(r^)  Is  a non-zero  common  denominator. 
?(*£)  has  almost  constant  rank  while  d(r^)  is  non-zero  almost  everywhere  and  hence 
can  effect  the  rank  of  P(rf)  only  on  an  algebraic  variety  (since  the  division  of 
a matrix  by  a nan-zero  scalar  does  not  effect  its  rank.)  As  such,  our  Jacobian 
matrix  has  almost  constant  rank  Implying  that 

18.  5(rf)  - n - rank  (rf)] 


is  also  almost  constant.  The  proof  of  che  Theorem  is  therefore  complete* 

Given  the  theorem,  we  may  now  define  a global  measure  of  solvability  for  the 

fault  diagnosis  equation,  <5 , as  the  generic  value  of  d(r^) . I.  a.  the  value  5 (rf) 
cakes  on  for  almost  all  r^.  This  proves  to  be  a natural  measure  of  solvability 
since  it  indicates  che  ambiguity  which  will  result  from  mi  attempt  to  solve  che 
fault  diagnosis  aquations  in  a neighborhood  of  almost  any  failures.  Of  course, 
one  requires  some  sort  of  equation  solving  algorithm*®’ **  to  locate  a neighbor- 
hood of  an  actual  failure.  The  5 parameter,  however,  represents  a bound  on  the 


» 
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performance  of  any  such  algorltha.  Finally,  we  noea  that  slnca  $ la  Independent 

of  r^,  the  solution  of  the  fault  diagnosis  equations.  It  can  be  co^ruted  at  the 

tlae  the  systea  aid  Its  test  algorltha  are  developed  by  evaluating  5 (r)  at  a 

randoaly  chosen  generic  point,  say  r . In  turn,  this  paraaetar  a ay  then  be 

o 

employed  as  an  aid  In  the  choice  of  test  frequencies  and  test  points. 


TEST  FREQUENCY  SELECTION 


Adopting  the  measure  of  solvability,  &,  formulated  In  the  proceeding  section, 
it  remains  to  develop  an  algorltha  for  choosing  a set  of  test  frequencies; 
s^,  s ^ . ...  s^;  which  maximize  the  solvability  of  the  fault  diagnosis  equations 
(i.e.  minimize  5).  To  this  end,  let  denote  the  minimum  value  achieved  by 

5 for  any  set  of  test  frequencies;  s^,  s^,  ....  s^;  k - 1,  2,  ....  . Since  the 
possible  values  for  5 are  restricted  to  the  finite  set;  5 - 0,  1,  ....  n;  such  a 
minimum  is  assured  to  exist. 

The  following  theorem  gives  an  explicit  formula  for  computing  6 . while  Its 

am 

proof  yields  an  algorltha  for  choosing  a set  of  test  points  which  achieve  5.  . 

am 

Since  the  purpose  of  this  theorem  Is  to  formulate  an  algorltha  for  choosing  test 
frequencies,  the  theorea  Is  expressed  In  terms  of 


19.  vec[S(s,r) ] - I^d-Zd.r)!^)"1]  vec[z(s,r)] 

and 

20.  dvsc[^s»r?  }a(  [(1  + lii(1.z(.tr)Lii)-42(atr)]Li2)et(L2i(1.z(Str)Lii)-l)} 

idvscCZ(s,r)]/dr  } 

viewed  as  rational  fractions  In  s rather  that  In  terms  of  the  function  F(r) 


which  is  formulated  In  terns  of  an  a-priori  choice  of  test  frequencies. 
Theorea  2;  Let  Z. (s,r);  1 ■ 1,  2,  ...,  q;  be  rational  In  s and  r.  Then 


6^w  • n - col-rank 


dvecfS(s.r) 


» 


I 


t 


t 


9 


♦ 


9 


9 


u 

whan  n la  tha  dimension  of  eha  parameter  vactor,  r,  and  "col- rank"  danotea 

cha  ganarlc  number  of  linearly  lndapandane  coluans  of  eha  rational  matrix  [dev[S(s,r)]/dri 

over  eha  flald  of  coaplax  numbers.  Honovar,  3m<n  la  achieved  by  almost  any 

cholca  of  n-d^^  distinct  coaplax  frequended. 

Proof:  For  eha  a aka  of  bnvlty,  va  will  prove  eha  ehaona  only  for  eha  spadal 

casa  vh an  S(s,r)  la  a scalar  eranafar  function  (allowing  us  eo  drop  eha  "vac" 

erans  formation)  though  aasanelally  eha  saaa  proof  goes  eh  rough  In  eha  general  case 

modulo  soma  notaclonal  cosp  laxities . ^ Also*  since  eha  rank  of  eha  Jacobian  matrix 

Is  almose  conseane  le  suffldas  to  fix  eha  parameter  vector,  r,  at  any  generic 

point,  say  r . This  then  reduces  [dvec[S(s,r) ]/dr]  eo  a rw  vector  of  raclonal 
o 

fractions 


21.  R(s)  - [R^s)  Rjts)  ...  Rn(s)] 


when 

22.  R^s)  - [dvae[S(s,ro)]/dr^] 

and  our  problem  n duces  to  eha  verifications  of  ehe  fact  Chat  the  number  of 
linearly  lndapandane  coluans  of  R(s)  ovar  eha  field  of  coup  lax  sealan  Is  equal 
eo  eha  maxima  possible  rank  of  eha  complex  matrix 


23. 


R(S1) 

R(s2) 


R(s  ) 


W 

W 


W 


w 

R2^S2^ 


W 


W 

VS2) 


W 


= col  (R(s^)) 


over  all  possible  choices  of  the  complex  frequencies;  s. , s , ...  s.  ; k = 1,  2, 

X Z.  K 


« 
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Now,  clearly  if  some  column  of  R(s),  say  the  nth,  is  dependent  on  the  re- 
maining columns,  then 

2».  V*>  * £ =jRj<S) 


for  all  s.  Then  by  applying  24.  individually  for  each  s^ 


25. 


col  (R  (s.)) 
n i 


n-1 

l c.col(R.(s, )) 
i=l  J J 1 


for  any  possible  number  or  choice  of  the  s^.  The  rank  of  the  matrix  of 
Equation  23  is  therefore  less  than  or  equal  to  the  number  of  linearly  inde- 
pendent columns  of  R(s)  over  the  field  of  complex  numbers. 

To  prove  that  equality  can  be  achieved  with  an  appropriate  choice  of 
n^min  com?^ex  test  frequencies,  s.,  we  invoke  our  assumption  that  5(s,r)  is 
a scalar  transfer  function.  Without  loss  of  generality,  we  may  assume  that 
R^(s)  through  R^(s)  are  the  linearly  independent  entries  in  R(s)  over  the  field 
of  complex  numbers  in  which  case  we  must  show  that  there  exists  complex  fre- 
quencies s , s^,  . ..,  s.^  (k  = q in  this  case)  which  make  the  first  q columns 
of  the  matrix  of  equation  23.  linearly  independent. 

If  q = 1,  R^(s)  is  not  identically  zero  (since  otherwise  it  would  be 
linearly  dependent)  and  hence  for  almost  all  s^,  R^( s^)  #0.  As  such,  the 
columns  in  this  trivial  one  by  one  matrix  are  linearly  independent.  With  this 
as  a starting  point,  we  will  use  an  inductive  argument  to  show  that  the  theorem 
holds  for  all  values  of  q.  We,  therefore,  assume  that  it  has  been  shown  that 
for  q = p there  exist  complex  frequencies;  s^,  s2,  ...  , s^;  such  that 

the  matrix  r—  _ 


26. 


R. (s.  ) R-(s, ) ...  R (s  ) 

XI  i X pi 


R1(s2)  R2(s2^ 


yv 


\<v 


R„(s  ) 
2 P 


R (s  ) 
P P 


r 
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has  linearly  independent  columns  and  we  desire  to  show  that  there  exists  an 
s such  that  the  matrix 

p+1 


— 

— 

27. 

w 

W 

• • e R ( S - ) 

P 1 

» 

Vi(s)  = 

Rl(s2) 

W 

• ••  R (s„) 

P 2 

Vi(y 

t 

w 

vy 

•••  Rp^sp^ 

wv 

R1(s) 

32(s) 

R (s) 

p 

Vl(s> 

has  linearly  independent  columns  for  s = so+^.  3y  virture  of  our  assumption  that  S(s,r) 


a s 


v.a_ar  -oth  and  -?fl(3)  are  square  and  we  may  test  for  linear  independence 


of  the  columns  of  Rp+^(s)  by  computing  its  determinent.  Expanding  27.  in  co- 
factors along  its  bottom  row,  we  obtain 

p+1 

29.  det(R  , (s) ) = l (-l)P+:!+iA  .R.(s) 


3*1 


P+i.  33 


Since  R^  has  linearly  independent  columns  Ap+^ip+j_  t 0,  hence,  the  coefficiencts 
in  the  summation  of  equation  28.  are  not  all  zero  and  thus  by  the  linear  inde- 
pendence of  the  R^(s)  the  summation  is  not  identically  zero.  As  such,  one  can 

choose  almost  any  s which  will  make  the  determinant  of  R . , (s  , , ) non-zero  thus 
p+l  —p+1  P+1 

assuring  the  R^+^  has  linearly  independent  columns  when  its  rows  are  evaluated 
at  the  complex  frequencies  s^,  s^,  ...»  s The  proof  of  the  theorem  is  thus 

complete . 

Note  that  the  proof  of  the  theorem  yields  a natural  sequential  algorithm  for 
choosing  test  frequencies.  Moreover,  for  the  scalar  case  we  have  shown  that  the 
number  of  required  test  frequencies  is  exactly  (equal  to  the  column  rank 

of  the  Jacobian  matrix).  In  the  general  case  where  S(s,r)  is  not  a scalar,  the 

S 

number  of  required  test  frequencies  is  less  than  or  equal  to  n-<5  . . " 

min 
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Although  the  theorem  implies  that  one  can  randomly  choose  almost  any 

n-5  . test  frequencies  to  maximize  the  solvability  of  the  fault  diagnosis 
min 

equations,  the  result  does  not  take  cognizence  of  numerical  considerations.  Al- 
though no  theory  yet  exists  for  choosing  test  points  with  numerical  considera- 
ations  in  mind,  it  has  been  our  experience  that  the  "well  posedness"  of  the  fault 
diagnosis  equations  is  quite  sensitive  to  the  choice  of  test  frequencies.5  In 
most  of  our  experiments,  we  have  worked  with  real  test  frequencies  to  eliminate  the 
necessity  of  working  in  the  complex  plane.  On  the  other  hand,  m is  most  easily 
measured  when  values  of  s^  on  the  jw  axis  are  employed  whereas  it  has  been 
suggested  that  test  frequencies  symetrically  spaced  around  a circle  in  the  com- 
plex plane  might  yield  numerically  "well  posed"  equations. 

Although  the  measure  of  solvability,  5,  for  the  fault  diagnosis  equations 
is  dependent  on  the  choice  of  test  frequencies,  as  well  as  the  properties  of  the 

unit  under  test,  5 . is  determined  entirely  by  the  UUT;  its  comoonents, 
mm 

connections  and  accessible  test  points;  and  is  completely  independent  of  the 
test  algorithm  employed.  As  such,  5^^  may  be  taken  as  a natural  measure  of 
testability^  for  the  UUT  which  characterizes  the  degree  to  which  the  fault 
analysis  equations  can  be  solved  given  an  optimal  choice  of  test  frequencies 
and  solution  algorithm.  Moreover,  <5^^  may  be  used  as  an  aid  for  the  optimal 
selection  of  test  points.  * ’ To  this  end  we  may  choose  a set  of  test  points, 
from  several  options,  so  as  to  minimize  5^^.  Alternatively,  we  may  attribute 
a cost  to  each  input  and  output  test  point  and  then  choose  the  least  cost  com- 
bination of  test  points  which  yield  a specified  5 . . This  latter  process  re- 

min 

duces  to  a rather  straighfarvard  integer  programming  problem  and  is  thus 

4 5 

readily  automated.  * The  technique  is  illustrated  in  the  examples  of  the 
following  section. 
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EXAMPLES 


An  an  initial  illustration  of  th«  theory  consider  the  RC  coupled  amplifier 
with  inductive  load  shown  in  Figure  2.  Here  we  will  take  E^  to  be  the  only 


Figure  2:  RC  coupled  amplifier  with  inductive  load. 

test  input  but  we  will  initially  allow  Zq,  i^,  i^ , and  to  all  be  taken  as 

test  outputs  with  the  measure  of  testability,  S ■ , being  used  to  extract  a 

min 

reduced  set  of  test  outputs  from  these  options.  A component  connection  model 
for  this  circuit  is  given  by 


28. 
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Taking  our  vector  of  potentially  variable  component  parameters  to  be 
r = coKil,  L,  C,  R)  each  with  unity  nominal  value,  we  obtain  a nominal  trans- 
fer function  matrix 

30.  S(s,r)  = 

whereas  .our  Jacobain  matrix  evaluated  at  the  nominal  parameter  values  is  given 


Mow,  an  inspection  of  this  matrix  will  reveal  that  it  has  four  independent 

columns  over  the  field  of  complex  numbers  and  hence  if  all  four  possible  outputs 

are  used,  we  will  have  5^  = 0 implying  that  the  fault  diagnosis  equations  have 

locally  unique  solutions . On  the  other  hand , if  only  two  outputs , E and  i_ , are 

o C 

measured,  our  modified  Jacobian  matrix  will  reduce  to  the  first  and  third  rows 

of  the  matrix  shown  in  equation  31.  which  has  column  rank  3.  As  such,  if  we  only 

use  these  two  test  outputs,  we  obtain  S = 1 and  hence  the  solution  to  the 

nun 

fault  diagnosis  equations  will  be  characterized  by  a single  arbitrary  parameter . 
In  this  latter  case,  with  only  E and  i_  taken  as  test  outputs,  theorem  2 

O N# 

implies  that  dF /dr  will  have  rank  3 for  almost  any  choice  of  3 = n - S . test 

mm 

frequencies.  Ghocsing  s^  = 1,  sn  = 2,  and  s^  = 3,  we  obtain 


which  has  three  linearly  independent  columns  as  long  as  g(l)  i 0,  g(2)  i 3 and 


g( 3 ) / 0.  Indeed,  in  this  example,  any  two  of  the  three  frequencies  would  have 
sufficed  to  yield  three  linearly  independent  columns.  Mote,  for  scalar  transfer 
functions,  tr.eorem  2 implies  that  frequencies  are  actually  required  hut 

for  matrix  transfer  functions  fewer  frequencies  may  suffice. 


Of  course,  for  the  circuit  of  Figure  2,  we  have  a choice  of  some  15  combin- 
ations of  the  four  outputs  with  which  we  may  choose  to  work  for  the  diagnosis 

of  the  circuit.  The  resultant  5 . 's  for  the  various  combinations  of  outputs 

min 

are  £iven  in  table  l.5 

Finally,  with  the  aid  of  Table  1,  one  may  readily  develop  a test  point 

4 5 

selection  algorithm  for  our  circuit.  ’ For  instance,  if  we  desire  to  find  the 

smallest  set  of  outputs  which  yield  a 5 <_  1 an  inspection  of  the  table  will 

reveal  that  and  i^ , i^  and  ic , or  Eq  and  ic  are  the  optimal  choices . Of 

course,  if  one  attributes  a cost  to  the  various  outputs  (determined  by  the 
convenience  of  making  the  required  measurements ) , then  we  may  further  dis- 
tinguish between  these  three  possibilities.  For  instance,  if  voltage  measure- 
ments are  deemed  to  be  easier  than  current  measurements , the  combination  of 
i.  and  i„  may  be  excluded  with  the  decision  between  the  remaining  two  options 
being  dependent  on  whether  it  is  easier  to  measure  the  circuit's  input 

current  ( i„ ) or  its  load  current  ( iT  ) . 

v-  L 


i 
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Table  1:  Measure  of  testability  for  the 

circuit  of  Figure  2 using  various 
combinations  of  test  outputs. 

As  a second  example,  consider  the  one  stage  transistor  amplifier  shown 
in  Figure  3 with  the  AC  equivalent  circuit  of  Figure  4.  Since  it  is  clearly 
impossible  to  distinguish  between  failures  in  the  two  parallel  bias  resistors, 
Ra  and  R^,  these  two  resistors  have  been  combined  into  the  single  resistor, 
in  the  component  connection  model  of  equations  33.  and  34.  Taking  all  of  the 
component  parameters  as  potentiall  faulty,  r becomes  a 12  vector  composed  of 
C^,  r^,  ....  R^  and  as  before,  we  take  all  parameters  to  have  the  nominal  value 
of  unity. 


Figure  3.  One  Stage  Transistor 
Amplifier 


Figure  4,  Amplifier  Equivalent 
Circuit 


1 

„ vi 


Once  again  we  let  the  input  voltage  be  the  only  test  input  for  the  system  and 

we  take  V , I , V_  , and  I , to  be  possible  output  test  points.  The  resultant  6 j 

° C1  Ra  * . 

for  each  of  the  15  possible  combinations  of  these  output  terminals  is  tabulated  m 

in  Table  2. 5 


» 


From  the  table  it  is  apparent  that  no  single  test  output  suffices  to 


yield  a 5 . =0  (perfect  testability)  though  <5  . - 0 can  be  achieved  using 

nun  min 

two  test  outputs;  V and  I_  or  V and  I . 

o o e 

CONCLUSIONS 

Our  purpose  in  the  preceeding  has  been  to  formulate  an  analytic  theory  in 

support  of  the  intuitive  art  usually  associated  with  the  design  of  a test 

algorithm.  With  the  aid  of  the  techniques  developed  above,  we  believe  that  it 

will  be  possible  to  develop  an  automated  test  program  generation  (AT?G)  algor- 

4 5 

ithm  for  linear  systems.  ’ Indeed,  such  an  algorithm  could  be  readily  combined 
with  the  same  computer-aided  design  (CAD)  algorithm  used  in  the  system  design 

3 

process.  Given  the  component  connection  equations  such  an  algorithm  could  be 
employed  to  automatically  (or  interactively)  choose  test  points  and  test  fre- 
quencies and  generate  the  required  set  of  fault  diagnosis  equations.  These 
could  then  be  stored  cn  tape  and  supplied  to  the  automatic  test  equipment  (ATE) 
in  which  a faulty  system  would  be  tested  and  the  fault  diagnosis  equations  solved. 

Although  we  do  not  propose  to  discuss  the  actual  solution  of  the  fault 
diagnosis  equations  here,  it  should  be  pointed  out  that  by  assuming  that  relatively 
few  components  have  failed,  say  p<<  n,  it  is  possible  to  develop  specialized 
algorithms  for  the  solution  of  the  fault  diagnosis  equations  which  are  far  more 

7 j.1 

efficient  than  standard  equation  solvers  in  this  application.  ’ ’ These  are 

typically  derived  from  the  fault  simulation  algorithms  used  in  the  diagnosis  of  

digital  systems  and  may  naturally  be  classified  into  "simulation  before  test”  and 

"simulation  after  test"  algorithms.  Some  of  the  algorithms  are  discussed  in 

references  7,  9,  10  and  11. 

Finally,  we  note  that  as  formulated  above,  the  measure  of  testability,  6 . , 

min 

assumes  that  any  combination  component  failues  is  possible.  If,  however,  we 
assume  that  at  most  p<<  n components  fail  simultaneously,  the  ambiguity  in  the 


solution  of  the  fault  diagnosis  equations  may  actually  be  less  than  5 . . 

min 

For  instance,  in  the  example  of  Figure  3,  with  only  Vq  taken  as  an  output 
^min  = 3*  t*ie  fault  diagnosis  equations  can  be  solved  exactly  if  we  assume 
that  only  one  parameter  is  out  of  tolerance.'*'0  The  point,  here,  is  that  even 
though  the  solution  of  the  fault  diagnosis  equations  in  n-space  has  three 
arbitrary  parameters  when  the  solution  is  restricted  to  the  one  dimensional 
manifold  of  parameter  vectors  in  which  all  but  one  coordinant  are  nominal  it 
is  unique. 
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ABSTRACT 

An  algorithm  for  the- inversion  of  a continuously  parameterized  family  of 
sparse  matricies  is  formulated  in  terms  of  a differential  equation  characterizing 
the  evolution  of  the  sparse  L and  U factors  of  the  given  family  of  matrices. 
INTRODUCTION 

In  the  various  algorithms  used  for  the  analysis  and  design  of  large-scale 
circuits  and  systems,  the  problem  of  inverting  a continuously  parameterized 
family  of  sparse  matrices,  M(r),  is  often  encountered.1-5  In  frequency  domain 
analysis,  this  might  represent  a transfer  function  matrix  which  one  must  invert 
over  a specified  frequency  range  while  in  time  domain  analysis,  such  an  M(r) 
arises  in  the  form  of  the  Jacobian  matrix  for  the  system  equations1 * * *  which  is 
depenGent  on  some  potentially  variable  parameter,  r.  Typically,  one  inverts 
M(r)  at  a discrete  set  of  points;  r^ , i = 1,  2,  ....  n;  using  a sparse  matrix 
algorithm.'  Indeed,  the  more  efficient  algorithms  exploit  the  fact  that  the 
matrices  M(ri)  have  a common  sparsity  structure  allowing  much  of  the  compu- 
tational overhead  to  be  shared  by  the  n inversions.1 

An  alternative  to  repeated  inversion  is  the  continuations  algorithm5  wherein 
one  integrates  the  differential  equation 

1.  Z(r)  = -Z(r)(dfVdr)Z(r)  ; Z(0)  = M(0)-1 

I 

to  obtain  M(r)-1  3 Z(r).  While  the  integration  of  Equation  1 is  far  more 
efficient  than  repeated  matrix  inversion  for  small  matrices.  It  fails  to  take 

* advantage  of  the  sparseness  of  M(r),  thereby  rendering  the  technique  in- 
applicable in  a large-scale  systems  context.  The  purpose  of  the  present  note 

is  to  present  an  alternative  continuation  algorithm  which  combines  the  LU  factor- 

• ization  technique  of  sparse  matrix  inversion  with  Equation  1. 


I 
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LU  FACTOR  DYNAMICS 

Recall  the  standard  sparse  matrix  inversion  technique®  wherein  one  factors 
a matrix  into  the  form  M ■ LU  where  L is  lower  triangular  and  U is  upper  tri- 
angular with  ones  along  the  diagonal.  We  then  represent  the  inverse  matrix  in 
the  form  M"1  * The  key  to  the  technique  is  that  both  L and  U and  their 

inverses  will  be  sparse  if  M is  sparse  though,  in  general,  M'1  is  not  sparse. 

As  such,  one  may  store  and  manipulate  the  inverse  of  a sparse  matrix  via 
its  sparse  upper  and  lower  triangular  factors,  U_1  and  L_1 , even  though  the 
inverse  matrix,  itself,  is  non-sparse.  These  ideas  are  combined  with  the 
continuation  algorithm  concept  in  the  following  theorem.7  Here,  the  notation 
U[M]  is  used  to  denote  the  strictly  upper  triangular  matrix  obtained  from  M 
by  setting  all  of  the  entries  of  M on  or  below  the  diagonal  to  zero.  Similarly, 
7[M]  denotes  the  lower  triangular  matrix  obtained  from  M by  setting  all  of  the 
entries  above  the  diagonal  to  zero. 

THEOREM:  Let  X(r)  and  Y(r)  be  solutions  of  the  matrix  differential  equation 
X = -Xu[Y(dM/dr)X]  ; X(0)  = U(O)'1 

2. 

Y * -1[Y(dM/dr)X]Y  ; Y(0)  = L(O)"1 

Then,  X ( r ) 3 U(r)"1  and  Y(r)  * L(r)~^  where  M(r)~^  3 U(r)”^L(r)’^  is  the  LU 
factored  form  of  M(r)”^.  Note,  if  M(r)  and  dM/dr  are  sparse  then  every  matrix 
involved  in  the  integration  of  Equation  2 will  be  sparse.  Moreover,  the  inte- 
gration may  be  carried  out  with  the  aid  only  of  a matrix  multiplication  algorithm 
.plus  a simple  procedure  for  extracting  the  upper  and  lower  triangular  sub-matrices 
of  Y(dM/dr)X. 

Proof  of  the  Theorem:  First,  we  observe  that  if  Y(0)  is  lower  triangular,  then 
Y will  be  lower  triangular  and  so  will  Y(r)  for  all  r.  Similarly,  if  X(0)  is 
upper  triangular  with  ones  on  the  diagonal,  then  X,  being  the  product  of  an 
upper  triangular  and  strictly  upper  triangular  matrix,  will  be  strictly  upper 


triangular.  As  such,  X ( r ) will  be  upper  triangular  with  ones  on  the  diagonal  for 
all  r.  Thus,  X(r)  and  Y(r)  have  the  correct  form  and  it  remains  to  verify  the 
equality  M(r)"1  * X(r)Y(r).  Here, 

X ( r ) Y ( r)  * X (0) Y (0 ) + [ [X(q)Y(q) ]dq  =•  X(0)Y(0)  + [ [X(q)Y(q)  + X(q)Y(q)]dq 

J 0 J 0 

3 X(O)YCO)  + f (-X(q)u[Y(q)(dM/dq)X(q)]Y(q)  - X(q)1[Y(q)(dM/dq)X(q)]Y(q) Jdq 
J0 

* X(0)Y(0)  + f {-X(q)CY(q)(dM/dq)X(q)]Y(q)}dq 
J0 

= X(0)Y(0)  + [P  CX(q)Y(q)](dM/dq)CX(q)Y(q)]dq 
^0 

Differentiation  of  both  sides  of  Equation  3 with  respect  to  r then  results  in 


4.  [X(r)Y(r) ] * [X(r)Y(r)](dM/dr)[X(r)Y(r)] 


Finally,  a comparison  of  Equations  4 and  1 reveals  that  X(r)Y(r)  * 
both  X(r)Y(r)  and  M(r)"1  satisfy  the  same  differential  equation. 


M(r)’1  since 


EXAMPLE 


Consider  the  family  of  matrices 


5.  1 r 
M(r)  = 

_-l  1_ 

Here,  M(0)  is  lower  triangular  and  hence  has  the  trivial  LU-factorization 

6.  T 1 0*1  f 1 ol  [Tl  o“ 

M(0)  * - 3 L(0)U(0) 

-1  1J  [-1  10  1 


while 


0 1 
0 0 
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As  such,  we  have 


8. 


L(0) 


0 

1 


and 


9. 


U(O)'1 


Now,  upon  using  an  Euler  integration  formula  CZ(h)  = Z(0)  + HZ(0)]  we  may 
estimate  U(.l)"1  and  L(.l)"^  via  the  equalities 


U(.l)"1  : U(O)'1  + ( .1 )U(0)_1 


= U(O)'1  - (.l)U(O)"1  U[L(0)‘1M(0)U(0)"1] 


"l 

o" 

r°  ■' 

"l 

-1/10 

_0 

1_ 

[o  ° 

0 

1 

and 


l(.l)'1  = L(O)"1  + (.l)L(O)'1 


= L(O)'1  - 1 CL(0)“1M(0)U(0)"1 ]L(0)_1 


Multiplying  these  estimates  then  yields 


12. 


, , fa/ioo  -9/i oo" 

U( .1 )”'L( .1 )”'  = 

9/10  9/10 


which  compares  favorably  with  the  exact  inverse 
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The  error  here  is  due  to  the  approximation  inherent  in  the  numerical  integration  ' 
process  and  can  be  reduced  by  use  of  a more  accurate  integration  procedure.  Of 
course,  the  result  of  the  theorem  is  exact  and  the  computed  value  for  M(r)"1  will 
be  as  accurate  as  the  integration  process  employed. 
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There  has  been  considerable  work  dealing  with  the  topic  of  filter- 
ing for  problems  with  state  dependent  noise  [1-3].  As  well  as  being  of 


theoretical  interest,  the  topic  is  of  some  practical  importance  since 
many  systems  are  better  modeled  as  having  multiplicative  disturbances 
instead  of  additive.  One  example  occurs  in  the  momentum  exchange  method 
for  regulating  the  angular  procession  of  a rotating  space  craft  [4], 

There  is  a disturbance  which  depends  on  the  procession  rates.  Another 
example  occurs  in  the  design  of  phase  lock  loops  [2].  The  phase  insta- 
bility of  an  oscillator  described  in  rectangular  coordinates  appears  as 
white,  state  dependent  noise.  If  one  received  a signal  which  consisted 
of  a large  number  of  sinusoids  of  various  frequencies,  each  having  phase 
distortion,  then  one  would  have  to  build  a high  order  filter  to  recover 
the  signal  using  existing  ,-ethods. 

INTRODUCTION 

The  design  of  high  order  filters  is  often  problematic  from  the  view- 
point of  on-line  computation.  Therefore,  a number  of  researchers  have 
been  interested  in  designing  filters  of  reduced  order  [5-8].  It  often 
happens  that  one  is  only  interested  in  estimating  a lower  order  linear 
transformation  of  a state  vector,  and  it  seems  reasonable  to  attempt  to 
do  this  with  a lower  order  filter.  Design  of  the  filter  parameters  is  a 
fixed  configuration  optimization  problem  [8-10].  In  such  problems,  the 
structure  is  not  necessarily  optimal,  but  given  the  structural  constraints, 
the  parameters  are  selected  optimally.  It  is  interesting  to  note  that 
these  problems  often  have  non-unique  solutions  because  there  are  too  many 
free  parameters.  This  feature  can  be  used  to  obtain  filters  which  are 
easier  to  implement  than  well-known  techniques  such  as  Kalman  filtering, 


35 


36 


even  when  the  fixed  configuration  filter  is  of  full  order  [8].  In  some 
cases,  there  is  no  performance  loss  associated  with  the  alternative  linear 
filter,  [8],  [11]. 

In  this  paper  we  seek  to  extend  the  reduced  order  filtering  results 
developed  in  [8]  to  problems  with  state  dependent  noise.  The  problem 
is  similar  to  that  considered  in  [12],  however,  in  [12]  a discrete  system 
model  was  considered,  and  only  a single  stage/optimization  was  performed. 
Here  a continuous  time  problem  is  considered,  and  the  matrix  minimum 
principle  [13]  is  used  to  obtain  a solution.  Because  we  allow  a driving 
term  in  the  filter  to  remove  any  a-priori  bias,  it  turns  out  that  the 
problem  has  singular  arcs,  which  is  not  surprising  considering  previous 
works  [3],  [11]  in  the  area.  A very  nice  feature  of  the  work  is  that 
in  some  cases  only  linear  two-point  boundary  value  problems  are  obtained. 
These  can  be  solved  either  by  a direct  use  of  linear  systems  theory  or 
by  a Riccati  equation  technique.  Under  certain  circumstances  only  a 
single-point  boundary-val ue  problem  must  be  solved. 

PROBLEM  STATEMENT 

The  system  of  interest  is  assumed  to  be  modeled  by  the  I to  stochas- 
tic differential  equation 

dx(t)  = A(t)x(t)dt  + dw(t) 

+ E |'xi(t)-pi(t)]Gi(t)dv(t)  (1) 

i=l  L J 

where  x(t)  is  the  state  vector  of  d*minsion  n and  y(t)  is  the  mean  value 

of  the  state  vector.  The  disturbances  are  zero  mean  incremental  Wiener 

processes  with  covariances 

E{dw(t)dwT(t)}  = Q(t)dt 
Ejdv(t)dvT(t)j  = H(t)dt 


(2) 


1 


» 


t 


♦ 


» 


t 
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t 
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It  is  not  hard  to  show  [14]  that  the  mean  value  vector,  y satisfies 

du(t)  = A(t)u(t)dt  (3) 

The  initial  condition  for  (1)  is  random  with  known  mean  and  variance 

£ I <lt0)  | = »0  «> 

Vari*(to)l=  Po  <5> 

Equation  (4)  is  obviously  the  initial  condition  for  (3). 

The  observation  vector  is  also  corrupted  by  state  dependent  noise. 

d y (t)  = C(t)  x ( t ) dt+dv(t)  + £ ["xi(t)-yi(t)l  Mi(t)dv  (6) 

i = l L J 

In  (6),  y(t)  is  the  observation  vector  of  dimension  m,  dv(t)  is  the 
additive  measurement  disturbance,  and  dv(t)  is  the  multiplicative  dis- 
durbance.  The  vector  v(t)  may  be  large,  and  seme  of  its  elements  affect 
the  dynamic  model  through  the  terms  , while  others  affect  the  obser- 
vational model  througn  the  terms  . The  additive  disturbance,  dv(t) 
is  a zero  mean  incremental  Wiener  process  with  covariance 

E J dv(t)dv^(t)  | = R(t)dt  (7) 

The  terms  w(t),  v(t),  v(t)  and  x(t0)  are  uncorrelated. 

Only  a linear  transformation  of  x(t)  is  to  be  estimated,  i.e.,  it 
is  desired  to  estimate 

z(t)  = N(t)x(t)  (8) 

where  z ( t ) is  a vector  of  dimension  l4=n. 

The  estimate  of  z(t),  which  we  call  z(t)  is  constrained  to  be  obtained 
by  the  filter  equation 

dz(t)  = |V(t)z(t)  + g(t)J  dt  + K(t)dy  (9) 

The  vector  g(t)  and  the  initial  condition,  z ( tQ ) are  to  be  selected 
so  that 

E { e ( t ) \ = 0 Vt  =[t0,  tfj  (10) 
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where  e(t)  is  the  error  vector 


e(t)  = z(t)  - z ( t ) 


(11) 


The  matrices  F(t)  and  K(t)  are  then  to  be  selected  so  that  a quadratic 
performance  measure  t^. 

J = Ej  / eT(t)  Qe( t)dt  + eT(tf)Se(tf ) } (12) 

to 

is  minimized.  The  weighting  matrix  S is  assumed  to  be  positive  definite 
symmetric.  The  weighting  matrix  Q may  be  positive  definite  or  zero  and 
is  critically  important  to  the  solution. 

GENERAL  SOLUTION 

In  order  to  proceed,  it  is  convenient  to  develop  an  equation  for 

r ] 

the  error.  From  the  Ito  differential  rule  il5j  , it  is  seen  that 

dz(t)  = N(t)dx(t)  + N(t)x(t)dt  (13) 

Using  (6),  (9),  and  (13)  it  is  seen  that  the  differential  equation  of 
the  error  is 

de(t)  = dz(t)  - dz(t) 

de  = [ (NA-FN-KC+fl)  x-gdt+Ndw-Kdv 

L r n - n 1 

+Fedt+  IN  21  x G.  - k£  x.M.I  dv 


or 


i = l 


i = l 


(14) 


In  (14)  we  have  introduced  the  notation,  x = x-u.  From  (14)  it  is  seen 
that 


provided  that 


d_F{e(t)}=  F(t)E  | e(t) } 
dt 


g(t)  = (NA-FM-KC+N)uU) 


(15) 

(16) 
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If  furthermore 

z(to)  = N(to)u(t0) 


it  is  clear  that 

E { e(tQ)  }=  0 


(17) 


(18) 


From  (15)  and  (13),  one  can  see  that  (10)  is  satisfied  so  that  (16) 
and  (17)  are  appropriate  selections.  If  g(t)  is  selected  according  to 
(16),  the  error  differential  equation  can  be  written  as 

de  = (NA-FN-KC+N)xdt  + Ndw-Kdv 

n _ 

+ Fedt  + £ x.  ( MG . - KM . )dv  (19) 


The  eouation  for  x is 


dx  = Axdt  + dw  + 


i = l 


x-i-idv 


(20) 


Clearly  x and  e are  both  zero  mean  processes. 

If  (19)  and  (20)  are  put  in  1 equation,  it  is  easy  to  see  how 
the  second  moment  matrix  defined  as 


P(t)  4 

Pxxft) 

Pxe(t) 

A 

E{x(t)xr(t)| 

E|  x(  t)eT(  t)} 

^ex^) 

pee ( 

E]e(t)xT(t)| 

E|e(t)eT(t)[ 

propagates.  This  is  useful  since  the  performance  measure  (12)  may 
be  written  as 


(21) 


tf 

J = tr{/  QPee(t)dt  + SPee(tf)f  (22) 

fco 

If  one  has  the  appropriate  constraint  equation,  the  optimal  selection 
of  F(t)  and  K(t)  may  thus  be  solved  with  deterministic  theory  using  the 


matrix  minimum  principle. 

Equations  (19)  and  (20)  may  be  written  as 


» 


- AP  + P y.A  + Q +yi 

Y Y AA  * 


(29) 


Pgg  = [na-fn-kc+n] 

+ FP-  + P F 
ee  ee 


i 

[na-fn-kc+n] 


p + p 
xe  ex 

+ NQNT  + KRKT  + K ^K1 


- N *KT  - K ^2TNt  + N 


(30) 


and 


pxe  * Apxe  * pxx  (NA-FN-KCtN)1  + P^F1 


+ qnt  + - v/ 


(31) 


In  (29),  (30),  and  (31),  the  terms  '{/  , \\r„,  and are  defined  as 

l C 3 


n 

V 

j=i 


3 G. 
xxij  1 


G. 1 

J 


(32) 


% = 


n 

E p. 

i=i 
j=i 


xxij  Gi 


n 

3 -E 

J i = l 
j=l 


P M. 

XX  • • • 

*A1J 


="JT 


- M T 

= Mj 


(33) 


(34) 


The  term  Pgx  is  simply  the  transpose  of  Pxg . Clearly  P ^ can  be  calculated 
independently,  and  can  thus  be  regarded  as  a known  quantity.  The  problem 
is  to  select  K and  F so  that  (22)  is  minimized  subject  to  the  constraints 
imposed  by  (30)  and  (31). 

The  Hamiltonian  for  this  problem  is  then 

H = tr  [ OP  + P aT  + P \J  + P .\T  / (35) 

\ ^ ee  ee  ee  xe  xe  ex  ex  I ' ' 
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where  Aee , Axe , and  Agx  are  Lagrange  multiplier  matrices  associated  with 
pee»  pxe»  and  pex  respectively.  The  constraint  equation  for  Pex  is  in- 
cuded  for  symmetry. 

The  optimal  solution  for  the  gain  K(t)  is  obtained  by  setting  the 
gradient  of  H with  respect  to  K equal  to  zero.  This  leads  to  the  ex- 
pression for  K. 

-1 

K ’Vl|\e  (PexcT  * "V  t Aex  Cxx0'  [R  +%]  <36> 

where  the  required  inverses  are  assumed  to  exist.  The  Lagrange  multiplier 
matrices  satisfy  the  equations 

•:‘ee  = = - {Q  + AeeF  + FTAee}  (37) 

' pee 

and 

Axe  ="!L  * ' !(NA-FN-KC+N)T  Aee  + ATAxe  +.\xeF}  (38) 

5pxe 

The  matrix  Aex  is  just  the  transpose  of  Axe-  initial  conditions  for 
(29),  (30),  and  (31)  are 

W 5 = po  (39) 

and 

Pxe  (toJ  - Po  N (to)T;  pee  <‘0)  * * <V  Po  fl  (to)T  (40> 

The  terminal  values  for  (37)  and  (38)  are  as  required  by  the  transversal ity 

condition  applied  at  the  terminal  time 


Aee  <V  - S 

(41) 

Axe  (t^.)  =0 

(42) 

and 
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Notice  that  Agg(t)  can  be  computed  separately  without  solving  the  rest 
of  the  problem  if  F is  known  beforehand.  However  at  this  point,  we  have 
not  yet  determined  how  F should  be  selected.  It  will  be  seen  that  this 
depends  in  a critical  way  on  the  nature  of  Q.  We  will  consider  two 
different  classes  of  problems. 

CASE  I. 

In  this  case,  we  assume  that  Q = 0.  The  meaning  of  this  is  that 
the  quality  of  the  estimation  algorithm  is  only  important  at  the  terminal 
time.  This  may  make  sense  for  rather  a large  class  of  problems.  The 
reason  that  this  case  is  of  particular  interest  is  that  the  selection 
of  F does  not  affect  the  Hamiltonian,  so  that  we  are  free  to  select  its 
value  based  on  other  considerations. 

Consider  that  part  of  the  Hamiltonian  which  depends  explicitly 

on  F. 

H*  = tr  |F0+eTFT}  (43) 

where 

©=  (Pee-NPxe)  Aee  + (pex-NPxx)  Axe  (4«) 

From  (39)  and  (40)  it  is  clear  that0(tQ)  = 0.  If  it  can  be  shown  that 
0(t)  = 0 for  all  t in  the  interval  of  interest,  then  a singular  arc 
exists.  The  Hamiltonian  is  independent  of  F.  In  this  case,  one  does 
not  need  to  specify  F to  stay  on  the  singular  arc.  Differentiating 
©gives 


(45) 
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The  bracketed  term  in  the  above  is  zero  whenever  K is  chosen  optimally, 
i.e.,  according  to  (36).  Hence 

0(t)  = F(t)©(t)  -e(t)F(t)  (46) 

and  (46)  implies  that0(t)  = 0 for  all  t>tQ  since  0(tQ)  = 0.  The 
selection  of  F is  thus  not  a performance  factor.  It  may  be  selected 
a-priori  so  thatAge(t)  can  be  precomputed.  It  may  be  selected  so  as 
to  achieve  some  other  objective  such  as  reduced  sensitivity,  computational 
convenience  or  to  minimize  some  alternative  performance  measure  specifi- 
cally involving  F. 

When  one  thinks  about  it,  the  singularity  with  respect  to  F is 
not  particularly  surprising.  Clearly  two  different  filters  can  even  pro- 
duce the  same  output  at  a particular  time,  given  the  same  input.  What  is 
interesting,  is  that  this  fact  is  generally  overlooked,  and  as  the  example 
problem  will  show,  that  an  alternative  filter  structure  can  be  relatively 
easily  implemented. 

CASE  II. 

In  this  case  the  weighting  matrix,  Q,  is  a positive  definite  sym- 
metric  matrix.  When  one  develops  an  expression  for0,  the  result  is 

0=  F0-  0F  +0Q  (47) 

instead  of  (46),  where 

n*  NP«.  - pee  (48) 


Thus  unless  fl  is  zero,  a singular  arc  does  not  exist. 


f 


It  is  easily  seen  thatfl(t)  does  not  equal  zero  unless  F is  sleeted 
appropriately.  From  the  initial  conditions,  0(tQ)  = 0.  Taking  the 
time  derivative  of  fl  we  get 
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0=  Ffi+QFT  + (NPxx-Pex)  (NA-FN-KC+N)1 
-K  [rKT-  v/n1  + ^3KT-CPxe] 


Examining  the  last  equations  we  see  that  if 


(NA-FN-KC+N)  = 0 


0 = Ffi  + QF 


This  follows  from  the  fact  that  when  (50)  holds,  A xe ( t ) is  zero  for 

all  t in  the  interval.  Consequently  the  expression  for  the  gain  becomes 

-1 

K = [PexCT  + [r  + ^3]  (52) 


and  (52)  is  sufficient  to  have  the  last  term  in  (49)  be  zero.  In  view 
of  (51)  and  the  fact  that  fl(t0)  is  zero,  it  is  clear  that  Q(t)  is 
zero  for  all  te Jt0,  tfJ  provided  that  (50)  holds  and  that  the  gain  is 
selected  optimally. 

When  fl(t)  is  zero,  it  may  be  seen  that  the  orthogonality  require- 
ment is  met  in  a reduced  state  space,  i.e. 

n(t)  = N(t)Pxe(t)-Pee(t)  = E { [z(t)-e(t)]eT(t)f  = 0 (53) 

Since  z = z-e,  (53)  may  be  written  as 

E { z (t)  eT(t)  } = 0 (54) 
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so  that  what  we  have  required  for  singularity  is  that  the  error  and 
the  estimate  be  orthogonal. 

When  N is  the  identity  matrix  and  there  is  no  state  dependent 
noise,  the  result  is  the  Kalman  filter,  with  the  requirement  (50)  that 

F(t)  = A(t)-K(t)C(t)  (55) 

which  of  course  means  that  the  filter  is  of  full  order.  When  the  filter 
is  of  reduced  order,  and  N is  constant,  what  we  have  is  the  observer 
constraint  equation  [16] 

NA-FN-KC  =0  (56) 

In  general,  when  Q>0,  (50)  is  a necessary  condition  for  a singular  arc. 
Clearly  it  is  not  always  possible  to  select  F to  satisfy  (50).  In  such 


(59) 
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SPECIFIC  F SOLUTIONS 


In  the  preceding  section  we  have  shown  that  when  one  is  only 
interested  in  estimation  at  a particular  time,  the  selection  of  F may 
be  based  on  considerations  other  than  optimality,  so  that  one  may 
pick  it  prior  to  optimization.  Furthermore,  when  Q>0,  it  may  not 
be  possible  to  find  an  F which  results  in  a singular  arc.  In 
that  case  one  may  opt  to  select  F prior  to  optimization.  In  this 
section,  we  will  see  that  when  F is  selected  a priori,  the  two  point 
boundary  value  problem  which  must  be  solved  for. the  selection  of  K 
is  linear,  and  hence  relatively  easy  to  solve. 

Consider  substituting  the  gain  expression  (36)  in  (31)  and  (38). 
The  resulting  expressions  are 


and 


Pxe  = APxe  + Pxx(?1+NA-FN)T  + PxeFT  + QNT  + NT 

[(Pxx  cTtVT  Axe  '« 

. HT  * C pJ 


(60) 


Axe  = -(NA-FN+fl)TAee  - AT  Axe  - AxeF  + CT(R+  M^)’1 

•_<pxx  cT  ♦ VT  * <cpxe  * ’5nT>  'ee]  <6l> 

When  F is  known  a priori,  both  Pxx  and  Age  are  known  in  the  sense  that 
they  may  be  precomputed.  The  above  equations  are  then  seen  to  give  a 
linear  TPBVP  in  the  matrices  Pxe  andAxe.  The  solution  may  be  obtained 
in  a straight  forward  manner  using  linear  systems  theory,  or  alternatively 
by  assuming  that  the  elements  of  ,\xe  are  linearly  related  to  those  of 
Pxe , and  obtaining  a solution  involving  a Riccati  equation.  The  values 
obtained  for  Pxe  and  Axe  may  then  be  used  in  the  gain  expression  (36). 
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We  cannot  overemphasize  the  importance  of  the  fact  that  our  result 
is  a linear  TPB-VP,  since  it  is  reasonable  to  expect  to  solve  a linear 
matrix  TPBVP.  Often  a nonlinear  matrix  TPBVP  is  so  difficult  to  solve, 
that  the  utility  of  the  result  is  questionable.  We  shall  explain  pro- 
cedures for  solving  a linear  TPBVP  by  looking  at  a particularly  easy 
case  in  which  Age  is  a scalar  times  the  identity  matrix.  This  results 
when  both  F and  Q are  scalars  times  the  identity  matrix.  When  this  is 
true,  (60)  may  be  written  as 


pxe  = ^nPxe  + *-]_2  A*e  + 91  (6^) 

where 

Ln  • A + F - L*  C (63) 

L12  =-L*  (CPxx  A-1ee  (64) 

0X  « Pxx  (NA-FN+N)7  + QNT  +'V1  NT  - L*  V V (65) 

and  where 

L*  = (pxxcT  + V (R  + '^3 ) ' 1 (66) 


Equation  (62)  is  of  the  form 

Axe  = L21  Pxe  + L22  Axe  + °2  (67) 

where 

l2l  * CT  (R+  'ty'1  C Aee  (68) 

L22  * 'aT"F  + cT*-*T  (69) 

02  = -(NA-FN+N)1  Aee  + CT(R+  ^3)-1  ^ee  (70) 


and 
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Let  L be  the  matrix 


II 

_J 

Lu 

r 

CM 

—i 

. L21 

l22 

and  <t>be  the  associated  state  transition  matrix  which  can  similarly  be 
partitioned 


» 

(t>  = 

12 

-^21 

^22- 

Then  the  solution  to  (62)  is 

Pxe(t)  = *U  (t,  t0)Pxe(tQ)  +$12  (t,t0)  Axe(t0) 

+ J ^4^j(t,x)  D^(t)  (^»T)  D2(OjdT 

to 
and 

Axe(t)  = *21^ ’^O^xe^o^  + $22^’^0^  Axe^o^ 

+ J f$2^(t,T)  D^x)  + <{’22  ( t , T ) D2  ( X )J  dx 

fco  L 

Applying  (74)  at  time  t = tf  gives 

»xe(tf)  = 0 * ♦21(tt.t0)pxe(t0)  * •22(tt,t0)  *xe(t„) 

*ff  [*21^tf','MT'  * *22^f,T^2^]^T 
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We  can  solve  (75)  for  Axe  (t0)  and  substitute  the  results  in  (73)  and 
(74)  to  obtain  the  solution  for  all  te|t0,  t^J  . 

There  is  another  approach  which  is  probably  preferable  in  most  cases. 
We  assume  that  Axe  is  linearly  related  to  Pxe  by  the  relationship 


Axe 

(t)  = u(t)Pxe(t)  + B(t) 

(76) 

Differenatiating  (76), 

one  obtains  the  differential 

equation 

Axe 

" <e  *u[L„Pxe  + D!*Ll2  UPxe  * 

L12 

Bj  + B 

(77) 

A1 ternati 

vely, 

from  (67) 

V 

- 1 0 
u2 1'  *e 

+ L22  ,JPxe  + l22B  + °2 

(73) 

Equating 

(77)  and  (78), 

we  get  for  U 

U 

+ UL11  + 

UL12U  = l2^  i-22^ 

(79) 

and  for  B 

UL 12® 

+ 8 + UD^  = L 22B  + D2 

(30) 

Since  Axe 

(tf) 

= 0,  the 

terminal  conditions  for 

U and  for  B are 

U (tf)  = 0 

(81) 

B (tf)  = 0 

(82) 

The  Riccati  equation  (79)  and  equation  (80)  can  be  solved  backwards  in 
time  from  the  above  terminal  conditions.  The  optimal  gain  may  then  be 
expressed  as 

K =rPxeV+N,2  + Aee" 1 (UPxe^)T(PxxCT+y2)]  [r^3] 


(33) 
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and  Pxe  is  evaluated  as 

Pxe  = LnPxe  + L !2  [uPxe  + b]  + Dx  (84) 

The  matrices  Age,  Pxx,  U,  and  B must  be  evaluated  off  line,  however,  Pxe 
and  K can  be  evaluated  on  line  if  this  is  desired.  Most  likely  these 
would  also  be  evaluated  off  line  and  K stored  for  on  line  calculation 
of  z(t)  using  (9). 


t 


t 
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EXAMPLES 

The  first  example  we  shall  consider  is  of  the  category  discussed  in 
Case  II.  We  assume  that  A(t)  is  zero,  N(t)  = C(t),  and  that  there  exists 
an  F such  that 

FC  = C - KC  Vte[t0,  tfJ  (85) 

then  if  CC^  is  nonsingular 

F = [cCT  - KCCT]  (CCT)  _1  (86) 

The  filter  equation  is 

dz  = CCT  (CC1*)  zdt  + K jdy  - zdtj  (37) 

The  initial  condition  for  (87)  is 


J(t0)  • c(t0)  Uo 


(88) 


The  gain  is  of  the  form 

K(t)  = [px6TCT  + Cf2]  [R  + ?3] 


(39) 
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where  Pxg  is  the  solution  to 

Pxe  3 Pxe(CCT)'1[c£T-CCTKTJ  + (Q+f1)  CT-4'2  KT  (90) 

Alternatively  since  CPxe  = Pee>  it  may  be  desirable  to  evaluate  (89) 
as  -1 

K(t)  = [Pee  + C,2]  [R  + f3]  (91) 


where  Pge  is  the  solution  to 


ee 


Pee(t)  =[cCT(CCT)'1  - kJ  Pee+P 

+ CQCT  + krkt  + k?3kt  - cT  kt 

J 


CC* (CC* ) 


T>-1 


-k] 


- K’foV 


CflC 


(92) 


The  reason  (92)  is  appealing  is  that  Pee  has  fewer  elements  to  calculate 
then  Pxe . 

The  next  example  is  concerned  with  the  very  simple  problem  of 
estimating  a constant  having  zero  mean  and  variance  1 prior  to  observations. 
The  observation  is  of  the  form 


dy  = x dt  + dv  + Mx  dv 


(93) 


where  v and  v are  zero  mean  white  noise  with  covariance  parameter  1 
and  M is  constant.  We  are  interested  in  estimating  the  value  of  x at 
time,  T.  Hence 

J « E j e2  (T)}  (94) 


and  this  is  a problem  of  the  category  referred  to  as  Case  1. 
putational  convenience  we  select  F = 0.  The  TPBVP  then  is 


For  com- 


xe 


xe 


Y ^ Pxe  + ‘\e  ^ 

= y (pxe  + O 


xe‘ 


(95) 

(96) 


» 
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even  though  it  is  obviously  a full  order  filter.  The  authors  feel  that 
the  nonunique  property  of  optimal  linear  filters  for  certain  cases  is 
a feature  which  one  should  take  advantage  of. 

REMARKS  AND  CONCLUSIONS 

We  have  extended  the  results  of  [8]  to  problems  having  state  de- 
pendent noise  in  the  observation  and  dynamical  equation.  Control 
theoretic  methods  have  been  used  to  solve  the  problem,  and  optimal 
solutions  have  been  shown  to  correspond  to  singular  arcs.  Different 
solutions  result  when  there  is  an  integral  performance  measure  than 
when  only  estimation  at  the  terminal  time  is  important.  In  some  cases, 
we  have  seen  that  it  makes  sense  to  select  the  filter  matrix  ahead  of 
time  and  then  optimize  the  gain.  The  computational  algorithms  associated 
with  such  prior  selection  are  particularly  convenient.  There  are  no 
terribly  difficult  TPBVP's  in  this  approach  and  that  is  why  the  authors 
feel  that  it  is  practical  and  useful,  both  for  full  order  and  reduced 
order  filters.  The  amount  of  off  line  calculation  necessary  to  simplify 
on  line  filtering  appears  to  be  quite  realistic. 
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Consider  the  nonlinear  system 

n-1 

x(t)  = f(x(t))  + l u. (t) g . (x (t) ) , x (0)  = x_6  M, 

i=l  1 1 0 

where  M is  a connected  C°°  real  n-dimensional  manifold, 

f ,g^,  • • • are  linearly  independent  vector  fields  on  M, 

and  u^,...,u  ^ are  real-valued  controls.  Since  the  integral 

manifolds  of  g^,...,gn_^,  if  any,  are  real  (n-1) -dimensional 

submanifolds  of  M,  such  a system  is  called  a hypersurface  sys- 
tem. Suppose  U is  the  largest  open  subset  of  M which  is  reach- 
able from  Xg  and  suppose  U / M.  It  is  shown  that  the  boundary 

of  U is  a C~  real (n-1) -dimensional  submanifold  N of  M and  N is 
an  integral  manifold  of  the  vector  fields  g^  , . . . , gn_]_  • More- 
over the  restriction  of  f to  N must  assign  vectors  pointing  in 
the  direction  of  U.  Such  a U is  called  the  region  of  reach- 
ability for  the  system  with  initial  point  Xg . Many  ideas  here 

parallel  those  used  in  several  complex  variables  for  the  studies 
of  regions  of  holomorphy  and  of  uniqueness  of  analytic  continua- 

• ® n 

tion  for  the  CR-functions  on  a C real  hypersurface  in  £ , n > 1. 


INTRODUCTION 

Let  M be  a connected  C*  real  n-dimensional  manifold 
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f,  g^»...,gn_^  be  C linearly  independent  vector  fields  on  M, 
and  be  real-valued  controls.  The  system 


n-1 

x(t)  = f (x(t)  ) + l u.  (t)g . (x(t) ) , x (0)  = x 6 M, 

i=l  1 1 0 


is  called  a hypersurface  system  since  the  number  of  g vector 
fields  is  n-1.  We  know  that  the  reachable  set  of  this  system 
contains  an  open  set  in  M (see  arguments  in  [8]) , and  we  de- 
note by  U the  largest  open  subset  of  M which  is  reachable  from 
Xq.  This  set  U is  called  the  region  of  reachability  from  x^ , 

and  if  U = M,  the  system  is  controllable  from  x^ . If  U / M 

then  we  prove  that  the  boundary  of  U is  a C°°  real  (n-1) -dimensional 
submanifold  N of  M and  N is  an  integral  manifold  of  the  vector 
fields  In  addition,  the  restriction  of  the  vector 

field  f to  N must  point  in  the  direction  of  U. 

This  article  is  arranged  in  the  following  way.  In  section 
2 we  give  definitions  and  relevant  examples.  Section  3 contains 
a local  theory  concerning  the  boundary  of  U under  the  assumption 

that  this  boundary  is  near  one  of  its  points.  In  section  4 
we  state  a theorem  from  [6]  concerning  a subbundle  of  the  tangent 
bundle  to  M and  allowing  us  to  remove  the  C1  restriction.  Then 
we  prove  our  main  result.  Theorem  4.2,  and  give  several  applica- 
tions. 
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DEFINITIONS  AND  EXAMPLES 

We  shall  use  the  classical  Frohenius  Theorem  and  Chow's 
Theorem  [2].  For  a statement  of  these  results  and  their  appli- 
cations to  control  theory  we  refer  the  reader  to  [1J . 

Of  interest  to  us  is  the  controllability  of  the  system 

n-1 

(2.1)  * (t)  = f (x  (t)  ) + l u.  (t)  g . (x(t) ) , x (0)  * X-  f M( 

i=l  1 1 0 

with  M,  f,  # . . . * and  u^ , . . . , un_^  as  in  the  introduction. 

We  let  T (M)  denote  the  tanaent  bundle  of  M with  fiber  T (M) 

x 

for  x € M. 

Recall  that  if  X is  a vector  vield  on  M,  then  a is  an 
intecra 1 curve  of  X if  a is  a C mapping  from  a closed  interval 
I CIR  into  M such  that 

= X(a(t))  for  all  t €1. 

^e^n3-tion  2.1  [8],  If  D is  a subset  of  T (M)  , then  an  integral 
curve  of  D is  a mapping  a from  a real  interval  (t,t'  ] into  M 
such  that  there  exist  t = tg  < t^  < . . . < t^  = t'  and  vector  fields 

X^,...,X^  in  D with  the  restriction  of  a to  [ ^ , t:^  ] being  an 

integral  curve  of  Xi , for  each  i = l,2,...,k. 

Definition  2.2.  Let  D be  a subset  of  T (M)  and  let  xQ € M.  A 
point  x 6 M is  D-reachable  from  xQ  if  there  is  an  integral  curve 
a of  D and  some  T > 0 in  the  interval  for  a such  that  a(0)  = Xg 
and  a(T)  =»  x.  A subset  S of  M is  D-reachable  from  Xg  if  every 
point  x € S is  reachable  from  xQ. 
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For  the  remainder  of  this  article  D is  the  subset  of  T(M) 

n-1 

given  by  the  vector  fields  f(x(t))  + £ u . ( t) g . (x  (t) ) . There- 
in 1 1 

fore  we  drop  the  D from  D-reachable. 

For  complete  vector  fields  (i.e.  vector  fields  which  can  be 

defined  for  all  -•  < t < ®)  it  can  be  seen  from  arguments  in  [8] 

(and  not  difficult  to  prove  for  the  special  case  of  hypersurface 

systems  like  (2.1)  even  if  the  vector  fields  may  not  be  complete) 

that  the  reachable  set  from  Xg  contains  a nonempty  open  subset  of  M, 
containing  Xg  in  its  closure. 

Definition  2.3.  The  largest  open  subset  U of  M which  is  reach- 
able from  Xq  is  called  the  region  of  reachability  from  Xg.  If 

U = M,  we  say  that  the  system  is  controllable  from  Xg. 

We  make  two  final  comments  before  introducing  some  illu- 
strative examples.  Since  we  are  considering  unbounded  controls 
(both  positive  and  negative) , no  generality  will  be  lost  in 
assuming  that  we  can  move  along  the  integral  curves  of  g^,...rg 
If  we  define  the  Lie  bracket  of  two  linearly  independent  vector 

fields  f and  g in  1R2  by  [g,f]  = ^f  - £^g,  then  we  find  that 

[g,f]  is  a linear  combination  of  f and  g.  This  means  that  there 
are  no  "new"  directions  in  which  a solution  to  the  system 

Jc(t)  = f (x  ( t ) ) + u(t)g  (x(t)  ) , x (0)  = xq€IR2 

can  move  (see  [1]). 

Example  1.  Consider  the  linear  system 


* 1 

’l  o' 

X1 

'o' 

S 

+ u ( t ) 

*2. 

0 1 

*2. 

i 

1 

■ f (x  (t) ) + u (t)  g 
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where  M is  the  open  right  half  plane  in  R . By  the  well  known 

2 

test  of  Kalman  [7]  this  system  is  not  controllable  for  all  of  R 
from  any  point  xQ  in  IR2  since  g ^ ® is  equal  to  ^ . Let 

Xq  = be  an  arbitrary  point  in  M.  The  integral  curve 

of  g through  (x^°,x2°)  is  given  by  the  vertical  line  x1  = x^0. 
Suppose  U is  the  open  set  in  M defined  by  { (x^ , x2)  | x^  > x^  } . 

Since  we  can  reach  any  point  on  x^  = x^,  our  only  hope  in  escaping 
from  the  set  U is  that  there  is  some  point  on  the  line  x^  = x^0 
at  which  f is  in  the  direction  of  the  complement  of  U.  However, 


f ( x ( t ) ) = on  this  line  and  x.,  >0,  implying  U contains  the 

2 I ^ 

reachable  set  from  xQ.  Letting  u(t)  = 0,  we  reach  every  point 

0 

x 

on  the  line  x_  * — g x^  to  the  right  of  (x..  ,x_°).  Using  infinite 

X1 

controls  and  the  interval  curves  of  g,  we  find  that  the  region 
of  reachability  of  our  system  from  Xg  is  the  set  U. 

Example  2.  Consider  the  linear  system 


0 l x 


1 O x 


+ u(t)  = f(x(t))  + u(t)g, 

1 


where  M is  the  open  upper  half  plane  in  <R  . From  Kalman  [7]  we 

2 

see  that  this  system  is  controllable  for  all  of  R from  any  point 


- *2  . 

Xg  6 R since 


f0  ll  To]  fll  , 

[l  oj  [1,-LoJ  » 


is  not  a linear  combination  of 


0 0 


Let  Xg  = (xx  ,x2  ) be  any  point  in  M.  Again,  the  integral 
curve  of  g through  (x1°,x2°)  is  given  by  the  vertical  line 


*2  = an<^  we  must  restrict  to  be  positive  for  points  on 

this  line  in  order  to  stay  in  M.  Let  U be  the  open  set  in  M 

given  by  { (x^x^  | x^  > x^  , x2>0}.  The  vectors  that  f assigns 

at  each  point  of  the  line  x^  - x^°  in  M are  in  the  direction  of 

U since  x2°>0.  This  is  also  true  for  each  vertical  line  in  M 

to  the  right  of  x^  = x^0.  Hence  the  integral  curve  given  by 

f (x  ( t) ) with  initial  value  (x.°,x-°)  must  move  forever  to  the 

right.  Those  integral  curves  of  g including  and  to  the  right 

of  x^  = x^0  intersect  the  integral  curves  of  f transver sally 

(defined  in  section  3)  at  each  of  its  points.  Thus  the  reach- 
able set  from  Xg  in  M is  U.  Since  we  cannot  move  outside  of 

M,  U is  the  region  of  reachability  from  xQ. 

We  remark  that  we  are  assuming  f and  g are  linearly  in- 
dependent on  M in  our  theory.  Suppose  that  we  relax  this 

2 

temporarily  and  let  M * R in  our  Example  2.  Starting  at  any 
point  xQ  = (x1°,x2°)  in  R2 , we  can  always  move  vertically 
(both  up  and  down)  along  x^  3 xi^'  so  we  can  assume  x2^  > 0. 

An  argument  as  in  our  previous  discussion  shows  we  can  reach 
the  set  { (x^Xj)  |x^  > x^0}  since  we  no  longer  must  remain  in 

• I 

the  open  upper  half  plane.  We  take  a new  point  (x^  ,x2  ) 
*0  ' 

with  x^  = x^  and  x2  <0.  Since  the  vectors  f assigns  along 
the  line  x^  ■ x^  are  negative  in  the  x^  sense  when  (x1°,x2) 
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is  in  the  open  lower  half  plane,  we  can  also  reach  the  region 
{ (xi ' x2 ) I X]_  ^ xj_  »X2  ^ Moving  along  the  integral  curves 

of  g will  give  us  the  entire  plane,  as  predicted  by  the  Kalman 
theory. 

Examples  1 and  2 suggest  a possible  solution  to  our  problem 
for  2-dimensional  linear  systems.  If  u is  the  region  of  reach- 
ability from  Xq  in  M,  then  the  boundary  of  U in  M should  be  the 

integral  curves  of  g on  which  the  vectors  given  by  f point  in 
the  direction  of  U. 

Since  we  are  most  interested  in  nonlinear  systems,  we 
examine  the  following  bilinear  system. 

Example  3.  Consider  the  system 


= f(x(t))  + u(t)g(x(t)), 

where  M is  the  set  (R2-{(0,0)}.  Brockett  [1]  states  that  if 

xQ  = (x1°,x2°)  is  in  the  positive  quadrant,  then  the  region 

of  reachability  U from  x^  is  contained  in  this  quadrant. 

The  integral  curve  of  g through  a point  on  x^  = 0,x2  > 0 is 

the  line  (0,x2>  with  x2  > 0.  Moreover,  the  integral  curve 

of  g through  a point  on  x2  = 0,x^  > 0 is  the  line  (x^,0) 

with  x^  > 0.  These  together  form  the  boundary  of  the  first 

quadrant  in  M.  The  vector  field  f assigns  vectors  to  this 
boundary  which  point  toward  the  first  quadrant.  Thus  there 
is  no  hope  of  a solution  starting  in  this  quadrant  to  leave  it. 


oo 


We  could  give  many  more  examples  at  this  time,  but  they 
would  all  hint  at  the  same  conclusion.  In  the  system  (2.1), 
the  important  items  to  check  appear  to  be  the  integral  mani- 
folds of  g^,...,gn  if  any  exist,  and  the  direction  of  the 

vector  field  f on  these  integral  manifolds.  We  next  examine 
these  conditions  for  regions  of  reachability  with  boundaries. 

C1  BOUNDARIES 

Let  U be  the  region  of  reachability  of  the  hypersurface 
system  given  in  (2.1).  Let  x be  an  element  of  the  boundary  of 
U (denoted  by  3U) , and  assume  3U  is  in  some  open  neighborhood 
W of  x in  M.  As  just  mentioned  we  need  to  consider  the  direc- 
tions of  f on  W D3U  and  the  possibility  of  having  an  integral 
manifold  of  g^,...,g  ^ through  x.  Recall  that  a differentiable 

submanifold  S of  M is  an  integral  manifold  of  g^,...,gn_^  if 

Ty(S)  is  the  space  spanned  by  g^,...,gn_^  at  y for  each  y£S. 

For  a more  thorough  discussion  of  integral  manifolds  we 

2 

must  consider  the  Lie  bracket  which  we  defined  earlier  for  (R 

(a  very  special  case) . Let  g^  and  g^  be  different  vector 

3g  3g. 

fields  on  M and  define  [g.,g.]  = g.  - -r — Ug . . This  may  or 

J.J  oXj  d X X 

may  not  give  us  a "new"  direction  in  which  to  move  (see  [1]). 

Let  La  be  the  smallest  Lie  algebra  generated  by  taking  suc- 
cessive Lie  brackets  of  the  g^,...,gn_^  given  in  equation  (2.1). 

If  we  get  a vector  space  of  the  same  dimension  at  each  point  of 
M,  then  Lft  is  a fiber  bundle  with  constant  fiber  dimension  n-1 
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or  n.  Moreover,  L^  is  a fiber  subbundle  of  the  tangent  bundle 
to  M. 

The  following  definition  is  essential  to  our  work.  Let 
and  S2  be  submanifolds  of  M of  dimensions  k and  n-k 

respectively. 

Definition  3.1.  The  manifolds  and  S2  intersect  transver sally 

at  a point  y €S,P,S2  if  and  only  if  Ty(S1)  © Ty(S2)  = T (M)  . 

Here  ffi  denotes  the  direct  sum. 

We  now  prove  a result  under  the  assumption  that  locally 

our  open  set  has  a boundary. 

Theorem  3.2.  Let  0 be  an  open  set  in  M which  is  reachable  from 
x-  for  system  (2.1) , and  let  x be  an  arbitrary  point  in  30. 

Suppose  there  is  an  open  neighborhood  W of  x in  M such  that 
W f|30  is  a C1  real  (n-1)  -dimensional  submanifold  of  M.  If  any 
one  of  the  following  conditions  holds,  then  0 is  not  the  region 
of  reachability  from  x^ . 

i)  the  fiber  dimension  of  L^  at  x is  n, 

ii)  the  integral  curve  of  some  gi#  1 < i < n-1,  is  trans- 
versal to  30  at  x, 

iii)  f assigns  at  x a vector  pointing  in  the  direction  of 
the  complement  of  0. 

Proof.  If  i)  holds  then  the  fiber  dimension  of  L^  at  all  points 

in  some  open  neighborhood  of  x in  M must  be  n (since  n is  max- 
imal). Thus  30  cannot  be  an  integrable  manifold  of  g ,...,g 

1 n-1 


near  x by  the  Frobenius  Theorem,  and  there  exist  a point  y € 30 


arbitrarily  close  to  x and  a g^,  1 < i < n-1,  such  that  the 
integral  curve  of  is  transversal  to  30  at  y.  Hence,  i)  re- 
duces to  ii)  . 

Next  we  assume  that  ii)  is  true.  If  the  integral  curve 
of  g^,  chosen  arbitrarily  from  the  set  g^,...,gn_^,  is  trans- 
versal to  30  at  x,  then  it  is  transversal  to  30  in  W H 30  / 

W being  an  open  neighborhood  of  x in  M (this  W may  be  a smaller 
set  than  our  original  W) . Following  the  integral  curves  of 
g^  that  start  in  WHO,  a reachable  set  from  Xg , and  continuing 

past  W D 30 , we  have  that  0 cannot  be  the  region  of  reachability 


:rcm  xQ. 


If  iii)  holds  at  x,  then  it  holds  for  all  points  in  W 030, 
and  the  argument  given  in  ii)  with  g^  replaced  by  f implies  the 

desired  result.  Q.S.D. 

It  is  interesting  to  note  that  condition  i)  does  not  depend 
on  W0  30  being  a manifold. 

We  seek  a minimum  number  of  necessary  conditions  that  an 
open  set  U c M be  the  region  of  reachability  from  Xg. 

Theorem  3.3.  Let  CT  be  the  region  of  reachability  from  x.  of 

the  system  (2.1).  Suppose  3u  is  a manifold  for  an  open 
neighborhhood  W of  oc£3u  in  M.  Then  W|H3U  is  an  integral  mani- 
fold of  g^, . . . ,g^_^  • and  the  vector  field  f assigns  to  W0  3U 

vectors  pointing  -in  the  direction  of  U. 

Proof . It  follows  from  part  ii)  of  the  preceding  theorem  that 
Wfl  3U  is  an  integrafl  manifold  of  g^,...,g  Hence  Wfl3U  is 


actually  a C*  submamifold  of  M.  Since  f,  g^r...,gn_^  form  a 
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linearly  independent  set  on  M and  WflBU  is  an  integral  manifold 
of  g^,...,gn_^,  part  iii)  implies  the  statement  concerning  f.  Q.E.D. 

We  shall  prove  in  the  next  section  that  the  hypothesis  3U 

is  near  x is  superfluous. 

THE  MAIN  RESULT 

The  following  theorem  was  proved  in  [6]  for  use  in  the  unique- 
ness of  analytic  continuation  problem  for  CR-distributions  on  CR- 
hypersur faces  in  <En,  n > 1.  The  statement  concerning  a C2 

boundary  can  be  relaxed  to  , or  we  can  simply  replace 
2 

by  C everywhere  in  the  preceding  section. 

Theorem  4.1  [6].  Let  M be  a C*  manifold  of  dimension  n, 

~~~~ 

and  let  H be  a subbundle  of  the  tangent  bundle  of  M with 

fiber  dimension  n-1.  Suppose  U -M  is  an  open  set  with  the 

2 

property  that  if  OCU  is  an  open  set  having  a C boundary, 
then  for  each  x 6 30  H3U  we  have  T^(30)  = (the  fiber  of 

H at  x) . Then  for  each  point  x € 3U,  there  is  a neighborhood 

00 

V of  x,  a real-valued  function  h € C (V)  with  nonzero  differ- 
ential for  all  point  in  V,  and  a closed  nowhere  dense  set 
EC  r such  that 

(1)  30  HV  = (x  € V|h(x)  €E}, 

(2)  for  each  J l €E,  = (x€v|h(x)  = 2.}  is  an  integral 

I » 

manifold  of  H;  i.e.  the  boundary  of  U is  foliated  by  integral 
manifolds  of  H. 

We  now  restate  Theorem  3.3  under  more  general  conditions. 
Theorem  4.2.  Let  U be  the  region  of  reachability  from  xQ  of 

the  system  (2.1).  Then  3U  * N is  a C integral  manifold  of 


L 


gl'***'gn-l  ^or  more  generally,  is  foliated  by  such  integral 

manifolds)  and  f assigns  vectors  on  N which  point  in  the  direc- 
tion of  U. 

Proof.  Let  H be  the  subbundle  of  T (M)  spanned  by  g^,...,gn_^. 

If  0 is  an  open  subset  of  U with  a C boundary,  then  an  appli- 
cation of  Theorem  3.2  and  Theorem  4.1  give  us  the  stated  con- 
clusion. Q.E.D. 

We  have  the  following  important  corollary,  the  proof  of 
which  is  obvious. 

Corollary  4.3.  Suppose  M is  connected  and  contains  no  integral 
manifolds  of  g^,  . . . ,gn_^  for  which  both  of  the  following  state- 
ments hold: 

a)  The  closure  of  the  integral  manifold  is  foliated  by 
integral  manifolds. 

b)  The  vectors  assigned  by  f on  this  integral  manifold 
always  point  in  the  same  direction  relative  to  the  integral 
manifold  (i.e.  if  this  manifold  divides  M into  two  components, 
the  vectors  must  point  toward  the  same  component) . Then  the 
system  (2.1)  is  controllable  from  any  x^ € M. 

4d 

Let  A denote  Hausdorff  measure  (see  [3])  in  dimension  d 
on  M.  Suppose  L is  the  set  of  points  on  which  the  Lie  algebra 
La  has  dimension  n— 1.  Then  L is  a closed  set  in  M,  and  the 

i 

Frobenius  Theorem  implies  that  L contains  the  integral  manifolds 
of  if  any  exist.  For  such  an  integral  manifold  we 

must  have  An  ^ (L)  > 0,  and  we  have  proved  our  next  result. 

r~»  — 1 

Theorem  4.4.  If  A (L)  = 0 and  M is  connected,  then  the  system 
(2.1)  is  controllable  from  any  xQ  CM. 


/ X 


Notice  that  if  M is  of  dimension  2,  we  always  have  integral 
curves  of  g for  the  system  x(t)  = f (x  ( t)  ) + u(t)g(x(t)), 
x(0)  = Xq.  Thus  Theorem  4.4  does  not  apply  in  this  case. 

We  state  two  theorems  from  [1]  and  indicate  in  a rather 
superficial  way  the  relation  of  these  theorems  to  this  present 
work.  We  restrict  our  attention  to  dimension  2 and  to  a 
hypersurface  system. 

Theorem  4.5  [ 1 ] . Suppose  f and  g are  vector  fields  on  a 

oo 

C real  2-dimensional  manifold  M.  Suppose  that  {f,g}  meet  the 
conditions  of  Chow's  Theorem  for  c”  vector  fields,  and  suppose 
that  for  each  initial  condition  xQ  the  solution  of  *(t)  = f(x(t)) 

is  periodic  with  a least  period  T(x^) . Then  the  reachable  set 

from  xQ  of  the  system  x(t)  = f(x(t))  + u(t)g(x(t))  is  the  set 

given  by  Chow's  Theorem. 

We  start  at  Xg  g M and  take  the  integral  curve  of  g through 

xQ.  Suppose  this  curve  divides  M into  two  connected  components 

and  M . If  the  solution  of  x(t)  = f(x(t))  starts  at  xQ  in 

the  direction  of  M+,  then  since  the  solution  is  periodic,  there 
is  some  point  on  the  integral  curve  of  g through  xQ  at  which 

the  vector  of  f is  in  the  M direction.  Of  course,  this  is  in 
keeping  with  Theorem  4.2. 

In  [5]  is  proved  a very  nice  generalization  of  the  following 
result,  which  we  state  in  dimension  2. 


r 
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Theorem  4.6  [5] . Consider  the  system 

x ( t ) = f ( x ( t ) ) + u(t)g(x(t) ) , x ( 0)  - Xq  / 

for  a C°°  real  2-dimensional  manifold  M.  Suppose  [f,g]  = ag 
on  M,  where  a is  a C*  function  on  M.  Then  the  reachable  set 
from  xQ  is  obtained  by  taking  the  integral  curve  a of  f 

through  xQ  (in  the  positive  time  sense)  and  then  all  integral 
curves  of  g interesecting  a. 

This  2-dimensional  version  can  be  seen  in  light  of  the 
following  result  found  in  [4],  The  one-parameter  group  of 
transformations  generated  by  f permutes  the  integral  curves 
of  g with  a change  of  parametrization  if  [f,g]  = ag  for  some 
C°  function  a on  M.  Interpreted  freely,  once  an  integral 
curve  of  f passes  through  an  integral  curve  of  g it  can  never 
return.  This  seems  to  be  in  agreement  with  Theorem  4.2. 

An  obvious  question  to  ask  is  if  the  necessary  conditions 
of  Theorem  4.2  are  also  sufficient. 

Theorem  4.7.  Let  x^tM  and  suppose  U is  the  smallest  open  sub- 
set of  M with  Xg£  U satisfying  3U  = N is  an  integral  manifold 
of  g^,...,gn_^  and  f assigns  vectors  to  N in  the  direction  of 

U.  Then  U is  the  region  of  reachability  from  xQ  for  the  system 

(2.1)  . 

) 

Proof . We  know  that  we  can  reach  an  open  set  and  by  the  theory 
developed  in  this  paper  we  have  that  we  can  reach  U.  The  im- 
portant fact  to  remember  is  that  to  leave  U we  must  break  through 
N near  some  point  x^£  N.  In  the  system 


73 


t 


n-1 

5c  C t)  = f(x,t)  + 1 u.  (t)g . (x(t)  ) at  the  point  x.  we  can  move 

i=l  1 1 1 

in  the  directions  f ,g. , . . . ,g„  . , -g.,...,-g  . , and 

i n- x i n- i 

n-1 

f + l u. (t)g.  for  the  appropriate  finite  u.'s.  Since  N is 
i-1  1 1 1 


an  integral  manifold  of  g^, 

with  i ^ j will  give  us  no 
from  x^.  Also,  since  f,g^. 


. ..,gn_^,  Lie  brackets  like  [g^,g^] 

"new"  directions  in  which  to  move 
...,gn_^  span  T(M),  the  brackets 


[ f , gi ] , i=l,...,n-l,  will  yield  vector  fields  which  are  linear 
combinations  of  f,g^,...,g  ^ (the  same  is  also  true  for  suc- 

cessive Lie  brackets) . The  only  linear  combinations  here  which 
we  can  use  at  x^  are  those  already  indicated  by  the  system.  Q.E.D. 


The  proof  of  Theorem  4.7  is  applicable  only  for  hypersurface 
systems  (or  certain  general  systems  that  behave  like  hyper sur- 
face systems) . We  shall  consider  results  as  in  this  paper  for 
general  systems  of  the  form 

m 

x(t)  = f(x(t))  + l u . (t) g . (x (t) ) elsewhere.  Such  a system 

i=l  1 1 

with  m < n-1  seems  much  more  difficult  to  handle  than  a hyper- 
surface system. 


» 


» 
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ABSTRACT 


A class  of  two-dimensional  recursive  digital  filters  called  symmetric 
half-plane  filters  is  discussed;  some  properties  of  these  filters  are  derived 
and  it  is  shown  that  in  certain  situations  these  properties  may  give  the 
symmetric  half-plane  filters  both  theoretical  and  practical  advantages  over 
previously  proposed  filters.  In  particular,  they  are  ideally  suited  to  highly 
parallel  processing. 

INTRODUCTION 

In  the  literature  on  2-dimensional  recursive  digital  filters,  two  main 
types  of  filter  have  been  studied;  these  are  the  quarter-plane  filter  (e.g.[l], 
[2])  and  tr.e  asymmetric  half-plane  filter  [3].  Basically,  the  two  correspond 
to  two  dif-'erent  concepts  of  causality.  The  general  stability  conditions 
for  a wide  class  of  filters  (including  symmetric  half-plane)  were  discussed  in 
[4];  unfortunately,  however,  those  filters  are  not  recursively  implementable  in 
general.  Here  we  will  consider  a class  of  filters  which  are  recursively  im- 
plementable, and  satisfy  the  same  stability  conditions  as  those  in  [4]. 

SYMMETRIC  HALF-PLANE  FILTERS 
I 

By  a symmetric  half-plane  filter  we  will  mean  a (causal,  recursive)  2-di- 
mensional digital  filter,  the  denominator  of  whose  transfer  function  is  of 
the  form 

A(Zr  Z2)  = 1 * J j Z”  0) 


2 


This  differs  from  the  filters  in  [4]  in  that  m goes  from  1,  rather  than  0, 
to  M;  i.e.,  this  filter  omits  all  of  the  row  m=0  except  for  the  constant 
term;  the  asymmetric  half-plane  filters  omit  half  of  this  row.  The  filter  (1) 
is  recursively  realizable,  since  the  computation  of  the  output  at  any  point 
depends  only  on  the  outputs  in  previously  computed  rows;  looked  at  from  another 
point  of  view,  each  row  of  output  depends  only  on  previous  rows  of  output.  This 
has  two  effects;  firstly,  it  focuses  attention  on  the  row  as  the  basic  element 
in  the  filter;  secondly,  it  implies  that  all  the  outputs  in  a given  row  may  be 

computed  in  parallel,  since  each  output  in  a row  depends  only  on  outputs  in 

previous  rows,  and  not  on  any  of  the  outputs  of  the  same  row.  This  is  the  main 

practical  advantage  of  this  class  of  filters  - it  would  be  of  significance,  how- 

ever, only  in  real-time  hardware  applications  of  2-dimensional  filtering,  and 
these  seem  to  be  few. 


SOME  PROPERTIES 


Using  the  methods  in  [4],  one  can  easily  derive  the  following: 

The  filter  (1)  (i.e.,  the  all-pole  filter  whose  denominator  is  A(Z^ .Zg ) ) is 
stable  if  A(Z1  .Zg)  J*  0 for  all  (Z^,  Z2)  such  that  | Z-j  [ * 1 and  |Z2|  < 1. 

We  note  that  this  set  is  the  same  as  that  for  the  symmetric  half-plane  filter 
in  [4];  it  is  smaller  than  that  for  the  asymmetric  half-plane  filters  [3].  It 
is  the  smallest  “instability  set"  (known  to  the  author)  of  any  recursively  im- 
plementable  class  of  filters. 

However,  there  is  a price  to  be  paid;  the  amplitude  response  of  the  filter 
is  restricted  as  follows: 

If  A(Z-j,Z2)  is  of  the  form  (1),  and  if  ^ is  the  transfer  function  of  a 


stable  filter,  then 
(2k 


log  | A(eJ 01 ,eJ02) | de2  = 0 


(2) 


independently  of  e*. 


/ y 


> 


Thus,  the  average  gain  along  any  line  of  length  2tt  parallel  to  the  e2-axis  is 

constant;  or  in  other  words  the  filter  cannot  have  variations  in  the  ej-direction 

in  overall  (average)  gain.  Equivalently,  if  the  cepstrum  of  |A(eJ01,  eJ02)|  is 

given  by  £ £ amn  Lz  , 

m--®  n=  -*• 

% 

then  a =0,  for  all  n. 
on 

This  follows  immediately  from  (2)  and  the  definition  of  the  cepstrum;  (2) 
will  be  proved  in  a forthcoming  paper. 

This  implies  that  in  order  to  realize  an  arbitrary  magnitude  function,  the 
filter  must  either  have  a (nonminimum-phase)  numerator  polynomial,  or  the  filter 
must  be  cascaded  with  a 1-dimensional  filter  in  . It  is  very  easy  to  calculate 
trie  ideal  amplitude  response  of  this  filter. 

DESIGN  AND  IMPLEMENTATION  CONSIDERATIONS 

It  is  conceptually  convenient  (and  in  a large  number  of  cases,  computation- 
ally efficient)  to  implement  the  Convolution  in  the  Z-j-direction  by  means  of  the 
Fourier  Transform.  (It  is  assumed  from  the  beginning  that  the  dimension  of 
the  array  to  be  filtered  is  a known  fixed  constant  in  the  Zj -direction,  i.e.,  each 
row  is  of  the  same  fixed  width).  From  this  point  of  view,  and  regarding  each  row 
as  a single  entity  described  by  its  1-dimensional  Z-transform,  the  coefficients 

a in  (1)  are  irrelevant;  what  matters  are  the  M functions 
mn 


Further,  the  stability  requirement  for  the  filter  is  equivalent  to  the  requirement 
that  for  each  fixed  9lt  the  filter  defined  by 


V 'z*>  - 


M 


I im(eJ9')Z?  * 1 


m=l 


is  1-dimensionally  stable.  Finally,  the  functions  am(eJ01)  do  not  have  to  be 
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analytic  or  meromorphic  functions;  this  is  seen  by  letting  N-«».  In  other  words, 
the  roots  of  am(eJ01)  can  vary  quite  arbitrarily  with  9X.  Thus,  we  can  design 
the  one-variable  filter  A0^(Z2)  by  any  one  of  the  usual  one-variable  design 
methods  we  choose  (yielding  a stable  filter)  for  each  the  result  will  be  a 
stable  two-variable  filter:  further,  if  our  one-variable  method  gives  poles  and 
zeros  explicitly,  we  have  the  same  for  our  two  variable  filter,  which  can  there- 
fore be  expressed  as  a cascade  of  filters  of  degree  1 in  Z2.  Finally,  if  one 
desires  a filter  of  finite  degree  in  Zj,  one  can  solve  the  following  approximation 
problem  (for  each  m,  1 <_  m ^M);  minimize  (over  bmn) 

II  £ bmneJn01  " bm(ej9l)ll  subject  t0 

n=-N 

N 

| y ejn9i|  < 1 for  all  9} , where 
n=-N  mn 

bm(eJ01)  denotes  the  m-th  pole  of  AQ ^ ( Z2 ) as  a function  of  elt  and  ||  |[  denotes 
some  error  norm. 

Hopefully  this  will  become  clearer  on  consideration  of  the  following-  example. 
EXAMPLE 

We  wish  to  design  a filter  with  second-order  Butterworth  response  in  Z2  to 
approximate  the  fan  filter  whose  passband  is  the  set  |e2|  < |9i|.  For  fixed  Sj, 
therefore,  the  filter  is  a 1-dimensional  filter  whose  passband  is  the  set  |e|  < | e x | . 
Using  the  bilinear  transform  technique,  we  find  the  second-order  continuous 
Butterworth  filter 

1 

1 + /2  s/iuc  + s2A>£ 


4 


transforms  into 


O X 


»C  0 + z2)* 

“c  ~Jl  “c  * 1 * <V  * J>  zi  ♦ («CJ  ♦ /*  »e  ♦ 1)  Z* 


while  its  (stable)  poles  transform  into 


(1  + "c(7^til,  }/{1  - “c(7TtJ)  1 

and 

11  + “c(7J~^)  l/l1  “ “C^7J“)  * 

for  <•>  >0. 

c — 


In  accordance  with  the  usual  frequency  warping,  we  take  uic  = | tan  9l/2[;  however, 
we  note  that  this  causes  stability  problems  at  9^0  and  9 j =-rr ; we  therefore  take 
a small  perturbation  ofuj.,  e.g.. 


'(cos  81/2)^  + £ 


The  filter  can  now  be  directly  implemented  by  multiplying  the  Fourier  transform 
of  the  previous  output  rows  by  the  appropriate  functions  according  to  (3)  and  (4) 
and  performing  the  recursion  from  row  to  row  directly. 

Alternatively,  if  a finite-order  (in  filter  is  desired,  we  must  solve 
the  problem:  Minimize  (over  bln) 


'•  -j.,  ■*  1 1 

subject  to  I * 1 for  3,1  •. 


and  similarly  for  the  other  root.  The  resulting  filter  may  then  be  implemented  in 
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cascade.  It  should  be  pointed  out  that  while  the  above  optimization  problem  is 
not  simple,  it  is  one-dimensional. 

Finally,  the  above  example  was  chosen  for  simplicity  and  convenience  rather 
than  realism.  Clearly  a Butterworth  filter  is  not  optimal  for  this  problem, 
especially  when  it  yields  a design  which  is  not  all-pole;  if  we  have  to  store 
input  rows  (as  well  as  output  rows)  we  may  as  well  use  them,  and  design  an 
elliptic  type  filter;  alternatively,  we  might  use  a filter  which  is  all-pole 
in  its  discrete  form. 
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ABSTRACT 


Many  physical  noise  processes  are  signal-dependent.  One  well 
known  example  is  film-grain  noise  (1-3).  In  this  note,  an  example 
of  the  application  of  optimal  estimators  for  images  in  signal-de- 
pendent film-grain  noise  is  presented. 

THE  MODEL 

A versatile  model  incorporating  both  signal-independent  addi- 
tive noise  and  signal-dependent  noise  is  utilized.  This  model  is 
given  in  Eq . ( 1) , 

r = s + kf (s ) + n_ , (1) 

where  r is  the  observed  photographic  density,  s is  the  original 
uncorrupted  image  density,  k is  the  scanning  constant,  f(s)  is  some 
function  of  s,  and  n^  and  n2  are  signal-independent  noise  processes. 
Thus,  the  middle  term  on  the  right-hand  side  of  Eq.  (1)  is  the  signal- 
dependent  noise  term. 

It  is  assumed  that  n^,  n2 t and  s are  mutually  statistically 

independent.  To  apply  the  model  to  film-grain  noise  problems, 

let  f(s)  = s^5,  where  p is  ususally  taken  to  be  1/2  or  1/3  (1-3). 

In  this  note,  we  let  p = 1/2  and  we  assume  n^  and  n2  are 

2 2 

zero  mean  Gaussian  random  variables,  with  variances  ct^  and  a^r  re- 
spectively. Further,  s is  assumed  to  be  a Gaussian  random  variable 

. . 2 
with  mean  u and  variance  a . 

s s 

THE  ESTIMATOR  STRUCTURES 

The  maximum  likelihood  (ML)  estimate  is  found  by  maximizing 
p(r/s)  over  s (3).  For  the  model  of  Eq.  (1),  the  estimate  is  found 
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to  be 


ML 


r2  + 


l2  ) k2o2  + \ k2a2/ J 


.2  2 

k 


.2  2' 
k 


(2) 


as  compared  to  the  simple  estimate 


SML  = r 


(3) 


which  results  when  the  signal-dependent  noise  term  of  Eq.  (1)  is 
omitted. 

The  maximum  a posteriori  probability  (MAP)  estimate  is  found 

by  maximizing  p(s/r)  over  s(3).  For  the  model  of  Eq.  (1)  and  the 

above  assunmtions,  the  estimate  s is  found  to  be  the  solution  of 

map 


.2  2 2 -.4  4 


+ p p-2* 

0 s s 

4 2 2 2 

[2u_-4k  o.OtP  . . - 

2 1 ~~*  * S * + 2»2> 

4 

+ jk2o2(o2-r2)  - 2a2r  - J = 3. 


(4) 


Again,  omission  of  the  signal-dependent  noise  term  in  Eq.  (1)  results 
in  a comparatively  simplified  estimate. 


3 MAP 


a 


2 

s 


2^ 

a +o 
s 


2 

2 


r + 


(5) 


Because  this  MAP  estimate  includes  prior  information  about  the 


image,  it  should  give  superior  performance.  In  fact,  under  the 
above  assumptions  it  can  be  shown  that  the  MAP  estimator  minimizes 
the  mean  square  estimation  error  (3). 

RESULTS 

Figure  1 is  the  original  image  of  an  archer.  Figure  2 is 

the  noisy  image  generated  digitally  according  to  the  model  of  Eq. 

(1).  The  image  in  Figure  3 is  the  estimate  found  by  the  solution 

2 

or  the  MAP  equation,  Eq.  (4) , with  u and  o taken  to  be  the 

s s 

sample  mean  and  variance  of  the  original  image. 

One  factor  severaly  affecting  estimator  performance  is  vio- 
lation of  the  assumption  that  the  image  statistics  are  Gaussian. 

For  a discussion  of  this,  see  the  paper  by  Froehlich  et . al . (3). 
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ABSTRACT 


Tn  the  following  wg  consider  a one-parameter  group  G acting 
on  a set  Q of  patterns.  We  show  that  under  extremely  mild  assump- 
tions the  problem  of  recognition  of  equivalence  of  two  patterns 
reduces  to  the  recognition  of  translated  functions  defined  on  the 
real  line.  Thus,  the  recognition  problem  is  reduced  to  one  of  a 
practical  nature.  In  fact,  the  well  known  properties  of  the 
Fourier  transform  with  respect  to  translation  can  then  be  applied 
in  most  instances  to  provide  a tractable  solution. 

FOURIER  TRANSFORMS 

Let  X denote  the  set  of  real  numbers  and  C denote  the  set  of 
complex  numbers.  We  follow  Bremerman  [1]  in  defining  the  Fourier 
transform  of  a function  f:  K - K (or  C)  as 


(a)  = / e-'  f (s)  as  . 

* — ao 


(1) 


Of  course,  the  integral  may  fail  to  exist,  even  for  v/ell-benaved 
functions  f.  Thus,  f may  be  a distribution  [1]  (i.e.,  a linear 

functional  on  a suitable  test  space).  In  any  event,  let  us  con- 
sider the  one-parameter  group  of  translations  acting  on  f : X ■*  K 
in  which  a point  t€  X corresponds  to  a translation 

f(s)  -*•  f (s-t)  . By  change  of  variable  in  equation  (1)  above 
we  have 


'A^(a)  = ejatf(a) 


(2) 


In  the  terminology  of  our  previous  reports  (6,7],  we  see  that  for 

A A 

each  a€C  the  mapping  E^:  f - f(a)  is  a relative  invariant  whose 
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modulus  is  the  exponential  nap  p^:  K -*•  { z C-  C | |z|=l)  gi^ven  by 

Da(t)  = e-jat.  This  is  a well-known  property  of  the  Fourier 
transform. 

Now  let  f^f^  K -*  K and  let  us  seek  t€  K such  that 

A A 

£2  ~ X^^i*  Provided  f^  and  f both  exist,  we  see  that  a necessary 

condition  is  that  for  all  a 6 c we  have 

f 2 fa)  = e^atf , (a)  (3) 

Determination  of  t now  becomes  a matter  of  practicality.  Observe 
that  the  condition  expressed  in  equation  (2)  is  sufficient  as 
well  since  it  implies  that  f^  = X a.e.  (almost  everwhere)  . 

Thus  we  see  that  the  values  of  the  Fourier  transform  give  a com- 
plete set  of  relative  invariants.  Now,  from  (2)  and  the  relation 

|eJ  | = 1 we  deduce  that  the  values  |f  (a)  | are  invariant.  Con- 
sequently, a necessary  condition  that  f2  = ^f^  for  some  t € K 

is  again  expressed  by  the  requirement  that  for  all  a € C we  have 

I f 2 (o)  | = | f L (a)  | . (4) 

In  the  following,  we  will  assume  that  a recognition  problem 
is  solved  when  it  is  reduced  to  the  translation  problem  for  func- 
tions defined  on  K. 
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ONE -PARAMETER  GROUPS 


By  a one-parameter  group  of  transformations  acting  in  Kn 
we  main  a family  G = {9?^}  of  coordinate  transformations  in  Kn 

which  depend  analytically  on  a real  parameter  t and  form  a 
group  with  respect  to  composition  (see  [4]  and  [5]).  Note  that 

each  v : Kn  ■*  Kn.  It  is  always  possible  to  choose  the  paramet- 
erization in  such  a wav  that  " *<F  = c?  ^ . Such  a svstem  of 

t s t + s 

coordinates  in  G are  called  canonical  coordinates . For  x e Kn 
and  t € K we  use  cp  ( x , t ) and  ^ (x)  interchangeably.  The  vector 

field  cf  the  group  is  defined  to  be 

Ux)  = lt=Q.  (5) 

Note  that  also  £ : Kn  -*■  Kn,  say  £ = ( ^ 2 , . . . , where  each 
Kn  ■*  K.  Finally,  the  infinitesimal  generator  of  G is  the 
differential  operator 

U = Z Z (x)|-  , (6) 

i-1  1 3xi 

where  x = (x, , . . . ,x  ) . 

l n 

Now  the  equations  for  canonical  coordinates  are  expressed 


cp  (fc(x,  t)  , s)  = cp  (x  , t+s) 

for  all  s , t £•  K and  x€Kn.  Application  of  ^ and  y^ls_o 


to  both  sides  of  equation  (7)  and  the  use  of  equations  (5)  and 
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(6)  yield  the  results 


U9  = 


So 

St 


and 


5 ('P)  . 


The  first  of  these  shows  that  O(x.t)  may  be  reconstructed  from 
the  infinitesimal  generator  as  the  unique  solution  to  the  partial 
differential  equation 


UO 


0 & 

a t 


cp(x,0) 


X , 


(8) 


while  the  second  shows  that  o(x,t)  may  be  reconstructed  from  the 
vector  field  t,  as  the  solution  of 


= 5 <«P> 

at 


<P(x,0)  = x, 


(9) 


which,  for  a given  x,  is  actually  an  ordinary  differential  equa- 
tion in  t.  Equation  (8)  is  typically  solved  [ ] by  obtaining 

solutions  to  the  system  of  ordinary  differential  equations 


dXi 

T^TxT 


dx. 


C2  (x) 


dx 

n 

5lV 


= dt . 


:io) 


Equation  (9)  may  be  expressed  in  integral  form  as 
CP( x , t ) = x + /q  E,  ( Cp(x,s)  ) ds. 


(11) 


The  latter  formulation  suggests  an  interesting  physical 
realization.  In  figure  (1)  we  see  a feedback  realization  of 
equation  (11) . The  output  of  an  integrator  is  input  as  argu- 
ment to  a generator  for  E,  values.  The  resulting  value  is  input 
to  the  integrator.  Starting  with  the  input  x to  the  integrator 
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at  tine  t = 0,  the  output  is  the  trajectory  cp(x,t),  t > 0. 
Replacing  £ by  -f,  gives  the  trajectory  backward  in  time, 
?(x,-t) , t > 0.  in  effect,  one  is  able  to  obtain  group  tra- 
jectories in  a straightforward  manner  in  analog  form. 

This  observation  merits  further  study  and  will  be  ac- 
tively pursued  in  the  future  research.  Comparison  with  the 
work  of  Brockett  [2,3]  amplifies  the  belief  that  further 
efforts  along  this  line  can  reduce  many  recognition  problems 


Figure  1.  Realization  of  equation  (11). 

to  problems  in  control  theory.  Initial  efforts,  not  complete 
at  this  time,  suggest  that  similar  methods  can  be  used  in  more 
general  n-parameter  groups  and  that  a control  may  be  introduced 
to  select  among  the  (infinitely)  many  one-parameter  subgroups. 
INVARIANTS  FOR  ONE  PARAMETER  GROUPS 

As  before,  let  G be  a one-parameter  group  of  transformations 

in  Kn.  It  is  shown  in  [5]  that  a map  h:  Kn  -*•  K is  G-invariant 
if  and  only  if  L'h  = 0.  Furthermore,  there  are  n-1  linearly  inde- 
pendent solutions  h , , h^ , . . . , h^_ ^ such  that  any  invariant  h is 

expressible  as  h = H (h, ,h_ , . . . , h .). 

l z n-l 


Let  ft  be  an  arbitrary  set  of  patterns  and  let  G be  a one- 
parameter  group  which  acts  on  ft.  Also,  let  R:  ft  -*•  V be  a mea- 
surement function  on  ft  with  values  in  some  real  vector  space. 
As  in  earlier  reports  we  obtain  a representation  of  w€ft  as  a 

function  wr : G ■»  V by  defining  wr(g)  = R(a  ^w)  . In  terms  of 
canonical  coordinates  in  G we  may  define  a new  representation 
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w:  K - V by  the  formula 

w(s)  = wr  (Z  ),  s^K.-  (14) 

s 

Now  for  t € K we  have 

c.  w(s)  = (cp.w)r(^  ) = wr(®  ~L-P  ) 
t t s t s 

= wr(cps_t)  = w(s-t)  = Xfcw(s)  . 

That  is,  we  have  successfully  represented  ft  as  a class  of 
functions  w:  K -►  V is  such  a way  that  the  action  of  G becomes 
translation.  The  question  of  equivalence  of  patterns  may  now 
be  resolved  by  the  use  of  the  Fourier  transform  methods  pre- 
viously discussed. 

CONCLUSIONS 

We  have  shown  that  the  recognition  of  equivalence  of  pat- 
terns under  a one-parameter  group  is  always  reducible  to  one- 
dimensional translation  of  functions  of  a real  variable.  The 
recognition  problem  is  thus  one  of  a practical  nature  which  may 
be  resolved  by  known  methocs,  such  as  the  use  of  Fourier  trans- 
forms . 

Further  work  will  be  done  in  the  attempt  to  generalize  the 
above  to  general  n-parameter  groups.  This  will  be  complicated, 
in  general,  by  the  lack  of  a transform  analogous  to  the  Fourier 
transform,  which  fails  to  be  readily  available  except  in  the 
Abelian  case. 

We  have  also  given  a means  of  analog  generation  cf  * 
tories  which  suggests  a control  theory  approach  to  rt> 
problems.  This  technique  will  be  extended  to  n-pa:  i- - 
and  the  relation  to  general  system  theory  will  be  -•  < 


REFERENCES 


[1]  Bremerman,  H. , Distributions , Complex  Variables  and  Fourier 

Transforms,  Addison-Wesley , Reading , Mass.,  1965. 

[2]  Brockett,  R.  W. , Non-linear  Systems  Theory,  Proceedings 

IEEE,  Vol.  64,  No.  1,  1976,  pp.  61-72. 

[3]  Brockett,  R.  W. , Lie  Theory  and  Control  Systems  Defined  on 

Spheres,  Siam.  J.  Appl.  Math.,  Vol.  25,  No.  2,  1973, 
pp.  213-225. 

[4]  Cohn,  P.  M. , Lie  Groups,  Cambridge  University  Press,  London, 

1957. 

[5]  Ince,  E.  L.  , Ordinary  Differential  Eauaticns,  Dover,  New 

York,  1956.  ‘ 

[6]  Newman,  T.  G. , Research  on  Pattern  Recognition,  in  Semi- 
annual Review  of  Research  under  the  Associate  Joint 
Services  Electronics  Program,  Contract  N0014-76-C-1136 , 
April,  1977,  pp.  96-103. 

Newman,  T.  G. , Pattern  Recognition,  in  Semiannual  Review 
of  Research  Under  the  Associate  Joint  Services  Elec- 
tronics Program,  Contract  N0014-76-C-1136 , October, 

1977,  pp.  143-160. 


[7] 


RESEARCH 


on 

BUILT-IN  TESTING 


R.  Saeks 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 


TEXAS  TECH  UNIVERSITY 


ABSTRACT 


A failure  prediction  algorithm  for  application  in  a periodic 
on-line  maintenence  system  operating  in  a Poisson  shock  environment 
is  described.  The  system  under  test  is  measured  at  periodic  main- 
tenence intervals  with  the  data  derived  therefrom,  being  used  to 
estimate  system  lifetime  and  determine  an  optimal  replacement 
time.  The  resultant  algorithm  is  simulated  and  compared  with 
various  fixed  replacement  schedules. 

INTRODUCTION 

Fault  analysis  processes,  have  been  and  will  continue  to  be  very 
significant  factors  in  the  safety  and  reliability  of  electrical 
systems.  This  is  especially  true  due  to  the  following  facts:  a 
rapid  advancement  in  the  complexity  and  size  of  modern  systems, 
increased  availability  and  capabilities  of  computers,  and  rapidly 
changing  technologies  in  integrated  circuit  fabrication.  Due  to 
this,  fault  analysis  has  become  much  more  than  an  academic  re- 
search topic.  Fault  analysis  is  applicable  in  an  industrial  en- 
vironment to  minimize  cost,  extend  the  lifetime  of  the  overall 
system,  control  maintenance  schedules,  and  effectively  plan  man- 
power needs. 

Although  considerable  effort  has  been  expended  during  the 
past  decade  to  develops  techniques  for  fault  detection  and  diag- 
nosis in  both  analog  and  digital  electronic  circuits*  little 
attention  has  been  given  to  the  possibility  of  formulating  algor- 
ithms for  fault  prediction.  To  accurately  predict  a fault,  a 
device  must  be  tested  at  periodic  maintenance  intervals.  If  the 
device  fails  or  does  not  operate  correctly,  it  is  replaced 
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immediately.  The  device  may  be  assumed  good  if  its  character- 
istics are  in  tolerence.  However,  if  the  characteristics  are  slightly 
off  nominal,  but  the  device  still  operates  correctly,  one  can 
attempt  to  predict  if  the  device  will  fail  before  the  next  scheduled 
maintenance  interval.  If  device  failure  is  predicted,  it  can  be 
replaced  before  failure  occurs  as  part  of  planned  preventative 
maintenance. 

With  the  advent  of  the  low-cost  microprocessor,  on-line  fault 

prediction  is  possible  and  practical.  A curve  fitting  algorithm 

for  on-line  fault  prediction  was  first  introduced  by  Saeks,  Liberty 
3 4 5. 

and  Tung  ' ' in  1975.  It  was  assumed  that  prior  life-time  statis- 
tics for  the  system  under  test  were  known.  Also,  performance  data 
of  the  system  at  each  maintenance  interval  were  collected.  The 
application  of  these  data  to  a second  order  polynomial  equation  re- 
sulted in  an  estimation  of  the  time  at  which  the  component  under 
test  would  exceed  tolerance  limits.  Based  on  a criterion  of  simul- 
taneously minimizing  on-line  failures  and  maximizing  component  life- 
time, a decision  as  to  whether  or  not  the  component  should  be  re- 
placed is  made  at  each  maintenance  interval. 

The  disadvantages  of  this  curve  fitting  algorithm  are:  the 
application  is  limited  to  failures  due  to  permanent  overstress,  the 
second  order  polynomial  is  too  simple  to  describe  the  performance 
of  the  component,  and  the  prior  lifetime  statistics  for  the  com- 
ponent are  often  not  available. 

Another  area  where  an  extensive  research  effort  is  being 

applied  is  shock  models  and  wear  processes.  Esary,  Marshall  and 
6 7 8 

Proschan  ' ' introduced  a shock  model  for  the  life  distribution 


of  a component  subjected  to  a sequence  of  shocks  randomly  occuring 
in  time  according  to  a homogeneous  Poisson  process.  They  also 
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considered  the  related  shock  models  in  which  each  shock  caused  a 
random  amount  of  damage  and  failure  occured  when  the  accumulated 
damage  exceeded  a specified  threshold.  This  failure  model  is 
well  known  in  modern  reliability  theory. 

Employing  the  Poisson-Shock  model,  another  curve  fitting  fault 
prediction  algorithm  which  will  overcome  the  disadvantages  of  the 

Saeks-Liberty-Tung  algorithm  will  be  discussed  in  the  present 
9 

paper.  In  the  following  section,  a model  for  the  failure  dynamics 
of  a system  component  parameter  is  formulated.  Here,  it  is  assumed 
that  the  faulure  is  due  to  the  component  being  subjected  to  a se- 
quence of  Poisson  distributed  shocks10,11  with  the  measurable  para- 
meter being  controlled  by  an  unknown  difference  equation  whose 
underlying  discrete  "component  time"  process  is  defined  by  the  number 
of  shock  to  which  the  component  has  been  subjected.  Since  both  the 
failure  dynamics  (i.e.  the  difference  equation)  and  the  relation- 
ship between  "component  time"  and  real  time  are  unknown,  our  fail- 
ure model  is  doubly  stochastic.  The  third  section  of  the  paper  is 
devoted  to  the  formulation  of  an  algorithm  for  estimating  the  com- 
ponent failure  dynamics  and  its  "lifetime",  defined  to  be  the  number 
of  shocks  required  to  cause  component  failure.  This  is  followed  by 
the  formulation  of  an  "optimal"  replacement  theory  wherein  the 
optimal  real  time  at  which  to  replace  a component  is  computed  in 
terms  of  its  estimated  "lifetime".  Finally,  the  results  of  a simu- 
lation of  the  algorithm  in  both  an  ideal  and  noisy  environment  are 
presented  and  compared  with  the  simulated  performance  for  several 
fixed  replacement  schedules. 

FAILURE  DYNAMICS 

The  performance  of  an  analog  device  subject  to  a series  of 
discrete  shocks  (switching  process,  improper  operation,  etc...)  may 
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drift  due  to  the  shock  damage.  Let  C(N)  represent  values  of  a 
particular  component  parameter,  where  the  "component  time",  N, 
denotes  the  number  of  shocks  the  component  has  received.  It  is 
assumed  that  the  drifting  parameters  can  be  described  by  a first 
order+  difference  equation  of  the  form: 

1.  C (N+l)  - C(N)  - aQ  - a^N  - a2N2  - ...  - ahNh  C(0)  =*  1 

Here,  the  coeficients  and  order  of  the  "forcing  polynomial"  are 
assumed  to  be  unknown  and  must  be  estimated  as  part  of  the  fault 
prediction  process.  A little  algebra  together  with  the  standard 
recursive  formula  for  solving  a difference  equations  will  reveal 
that 

N-l  h 

2.  C (N)  = 1 - l l a j1 

j=0  i=0  1 

Now,  if  the  tolerence  limit  for  the  component  parameter  is 
taken  to  be  C = 0 , we  may  define  the  lifetime  of  the  component  to 
be  the  smallest  integer,  N,  for  which  C(N)  0.  This  integer  which 
we  denote  by  L then  represents  the  number  of  shocks  necessary  to 
cause  the  component  to  fail. 

Consider  a simple  example  where  the  "forcing  polynomial"  is 
taken  to  be  of  the  first  order  with  positive  coeficients.  Then  1. 
reduces  to 


C (N+l)  - C (N)  -aQ  -ajN  ; C(0)  - 1 
From  equation  2.  C (N)  can  thus  be  expressed  as 


C(N)  » 1 - aQN 


N (N-l) 

2 al 


+The  concepts  described  herein  carry  over  without  modification  to 
the  case  where  the  failure  model  is  characterized  by  higher  order 
difference  equations.  The  first  order  model,  however,  suffices  to 
illustrate  the  theory  and  his  hence  used  throughout  the  present 
paper. 
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Then  the  lifetime  of  this  component  is  the  smallest  integer  satis- 
fying the  equation 


1 


a0L 


L(L-l) 

2 


a1  > 0 


» 


» 
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That  is,  L is  the  smalles 


L > 


(2a0  - ax)‘  * 


2a 


integer  such  that 

8al  ~ (2a0  ~ al} 
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Since  the  failure  model  of  equation  1.  is  dependent  on  "com- 
ponent time",  i.e.  the  number  of  shocks  the  component  has  received, 
rather  than  real  time,  it  remains  to  define  the  relationship  between 
"component  time"  and  real  time.  Following  common  practice  in  re- 
liability theory6 , we  assume  that  this  relationship  is  determined 
by  a Poisson  process.  Indeed,  this  is  the  unique  point  process 
which  has  the  scaling  properties  required  for  such  an  application.^ 
Here,  the  probability  of  N shocks  occuring  in  the  time  interval  t 
is : 

. 3.  PN(t)  - e"kt  - N =«  0,1,2,  ... 


^ Where  k is  a given  constant  representing  the  average  number  of 

shocks  per  unit  time.  Therefore,  (kt)  is  the  average  number  of 
shocks  in  the  time  interval  t. 

^ If  a component  with  lifetime  L is  subjected  to  Poisson  shock 

with  constant  k the  probability  that  it  will  fail  (i.e.  receive  at 
least  L shocks)  by  time  t,  is  then  given  by  the  formula9 

» 


xuu 


L-l 

F (t ) = l P(t) 

n=0 

l e"kt  (kt)n 
n=0  n! 

Thus,  even  though  the  lifetime  of  our  component  is  integer  valued, 
in  our  model  the  actual  failure  time  is  continuously  distributed 
since  the  time  at  which  the  component  receives  the  Lth  shock  is 
continuously  distributed. 

ESTIMATION  OF  FAILURE  DYNAMICS  AND  LIFETIME 

In  a periodic  maintenance  system,  the  performance  of  a component 

is  measured  at  each  maintenance  interval  nT.  That  is  to  say, 

(C^ , ^2 , •••/  Cg)  is  the  performance  data  taken  at  maintenance  times 

(T  , 2T , ...,  gT) . The  estimation  problem  can  be  stated  as: 

"Given  performance  data  (C^,  C2,  ...,  Cg) , T and  k,  estimate 

the  unknown  constants  (aQ,  a^  ...,  ah)  of  the  failure  dynamics." 

Since  it  is  assumed  that  the  system  is  subjected  to  Poisson  Shock 

with  constant  k,  the  expected  number  of  shocks  in  each  maintenance 

interval  is  kT.+  As  such,  if  we  assume  that  C is  the  value  of 

m 

the  component  parameter  at  N * mkT,  then  upon  substituting 
* C (mkT)  into  equation  2.  we  obtain 


t 


Although  not  theoretically  necessary,  we  assume  that  kT  is  an  integer. 
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Since  the  number  of  data  points,  g,  is  typically  much  greater 
than  the  order  of  the  polynomial  assumed  in  the  failure  model,  h,  it 
is  not  expected  that  equation  4.  admits  an  exact  solution.  Rather, 
we  attempt  to  solve  for  a coefficeint  vector,  A,  which  minimizes 
the  error  between  JA  and  Z.  In  particular,  if  one  adopts  a least 
squares  error  criterion  the  optimal  A is  given  by 


5. 


J_GZ 


where  J denotes  the  generalized  inverse  of  J.  Indeed,  if  as  is 

typically  the  case  J has  full  column  rank  than  J * (J  J)  J where 

"t"  denotes  matrix  transposition.  As  such,  we  take  the  A°  * 

col(a°,  a°,  ...,  a°)  as  our  estimate  of  the  coefficients  of  the 
o 1 h 

difference  equation  characterizing  the  failure  dynamics  of  our  drift- 
ing parameter , C , as  per  equation  1 . - - - - - 

To  estimate  the  failure  dynamics  of  a drifting  parameter,  the 
proper  choice  of  the  order  h is,  in  general,  quite  difficult  and 
depends  upon  physical  considerations  and  engineering  experience. 
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Once  h is  preselected,  however,  coefficients  to  best  approximate 
the  failure  dynamics  can  be  readily  computed  via  equation  5 . The 
accuracy  of  the  resultant  estimate,  however,  is  highly  dependent 
on  the  choice  of  the  order,  h,  and  on  the  number  of  measurements 
which  are  taken,  g.  To  find  a new  set  of  coefficients  for  a 
different  combination  of  h and  g,  the  entire  calculation  procedure 
is  typically  repeated  from  the  very  beginning  which  is  impractical 
in  the  on-line  maintenance  system.  Fortunately,  sequential  refine- 
ment schemes  for  obtaining  new  sets  of  coefficients  without  re- 

12 

peating  the  entire  calculation  can  be  developed.  ' As  such,  it 

is  possible  to  sequentially  update  onels  estimates  of  the  parameters; 

a , a^,  ...,  a.  ; as  additional  measurements  are  taken  and/or  to 

increase  the  order  of  the  model  for  the  failure  dynamics  without 

repetitious  matrix  inversion.  Our  algorithm  for  estimation  of  the 

failure  dynamic  underlying  the  measured  data  my  thus  be  readily 

implemented  on-line  with  the  computational  power  presently  available 

in  today's  microprocessors.  The  matrix  algebraic  details  of  the 

12  1 2 

required  sequential  refinement  schemes  are  straighforward  ' and 
readily  available  in  the  literature.  As  such,  they  will  not  be 
repeated  here. 

In  practice , given  g measurements  , C2 , . . . , Cg  taken  at 
maintenance  intervals  T,  2T,  3T,  ...,  gT,  one  sequentially  esti- 
mates the  coefficients  of  the  failure  dynamics;  aQ,  a^,  ...,  a^; 
increasing  h until  no  further  error  reduction  is  achieved.  The 
resultant  set  of  coefficients  is  then  used  in  equation  2.  to  de- 
termine the  component  lifetime,  L.  Upon  solving  the  equation. 


the  resultant  estimated  lifetime  is  found  to  be  the  smallest  in- 
teger, L,  such  that 

L-l  h 

l l M1  i i 

j=0  i-0  1 

Of  course,  if  the  measured  data  is  not  decaying  towards  zero, 
i.e.  the  component  isn't  failing,  this  inequality  will  have  no 

Q 

solution  in  which  case  we  take  L to  be  infinite. 

REPLACEMENT  THEORY 

Although  the  algorithm  outlined  in  the  preceeding  secion  yields 
an  "optimal"  estimate  of  the  number  of  shocks  required  to  cause 
failure  the  time  at  which  the  Lth  shock  takes  place  is  statistical 
in  nature  and  hence,  it  still  remains  to  determine  the  optimal  (in 
an  approprate  sense)  time  at  which  to  replace  the  component.  One 
such  criterion  is  formulated  in  the  following.  For  this  purpose, 
it  is  assumed  that  L has  been  computed  to  our  satisfaction  and  we 
desire  to  choose  a time,  Tr,  at  which  to  replace  the  component  as  a 
function  of  L.  Given  L and  Tr  we  denote  the  resultant  probability 


of  on-line  failure  (i.e.  failure  before  T^)  by  P 


r Pr 


1 - p. 


then  denotes  the  probability  that  the  component  is  replaced  at  time 
Tr  before  it  fails.  Similarly,  we  let  Tf  denote  the  expected  time 
to  failure  for  those  components  which  fail  on-line,  we  let  T de- 
note the  expected  time  to  failure  for  all  components  and  we  let  T* 
denote  the  expected  time  to  failure  for  the  components  if  they  were 
operated  to  failure  without  replacement  (i.e.  T*  = tJ _ ^ _)  . Finally, 
we  let  f^(t)  denote  the  probability  density  function  that  the  com- 
ponent receives  the  Lth  shock  at  time  t,  given  that  the  component 
fails  on-line,  whereas,  p^(t)  represents  the  density  function  of 
the  Poisson  distribution  with  parameter  (k)  and  E^ft)  represent  the 


112 


corresponding  distribution  function;  i.e. 

_ ...  _ (kt)1  -kt  „ , „ 

7.  p.(t)  - • i s i - 0|  1/  2/ 

L-l 

8.  E (t)  = l p.(t) 

i=0 

Q ^ * 

With  the  aid  of  some  elementary  calculus  anc*  T' 

as  well  as  their  derivatives  with  respect  to  T^.,  can  be  computed 
analytically.  As  such,  upon  defining  an  appropriate  cost  measure 
an  explicit  formula  for  determining  an  "optimal"  T given  L can  be 
derived.  We  begin  with  the  derivation  of  the  explicit  formula 
for  the  various  quantities  involved  in  our  replacement  theory. 

Since  a component  will  be  replaced  by  our  algorithm  if  an 
only  if  it  is  still  operating  at  time  T , i.e.  if  it  has  not  yet 
received  L shocks  at  time  , the  probability  of  replacement  is 
just  the  probability  of  receiving  less  than  L shocks  by  time  Tf. 

We  thus  have: 

(prop  1) 


(prop  2) 


Pr  = EL(Tr) 


Pr»  I 

i=Q 


L-1  -kT 

■ r—  e r 


L-l 

Jo  Pi<T^ 


EL(Tr) 


PC  * 1 - W 


ilJ 
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PL-l(t) 


(prop  4) 


fL(t)  " 


X/k  (1-EL(Tr)) 


To  derive  this  conditional  density  function  we  partition  the 
interval  (0,Tr)  into  N segments  of  length  A = Tr/N  and  we  compute 
the  probability  that  the  Lth  shock  takes  place  in  the  ith  time 
interval,  ((i-1) A,  iA].  Since,  this  can  be  caused  by  having  L-l 
shocks  before  (i-1) A and  at  least  one  shock  in  the  interval 
((i-1) A,  iA]  or  by  having  L-2  shocks  before  ( i— 1 ) A and  at  least 
two  shocks  in  the  interval  ( ( i— 1 ) A , iA],  etc.  the  probability  of 
failure  in  the  ith  interval  is  given  by 


l p . ( (i-1) A) [1  - E.  (A)  ] 
j = l J D 


I p ((i-UA)  l 
3=1  q=3  4 


(Ak) 


-Ak 


l p.  C l PL_S ( (i~l) A) ] (Ak)re"ak 


:=1 


s=l 


Taking  the  probability  density  function  at  a point  t in  the  interval 
((i-1) A,  iA]  to  be  limiting  value  of  this  quanity  divided  by  A as 
A goes  to  zero^  it  is  observed  that  the  terms  of  9 . containing 
powers  of  (Ak)  greater  them  1 go  to  zero  in  the  limit.  As  such, 
the  probability  density  function  for  the  Lth  shock  to  take  place 
at  time  t is  given  by 


10. 


limit  PL-1 ( Uk)e 


- k 


A-*0 


kpL-1 ( (i-1) A) 


1 


ilj 


> 


Finally,  since  we  are  only  interested  in  the  conditional  probabil- 
ity density  function  that  the  Lth  shock  will  take  place  at  time  t 
given  that  the  component  fails  on-line  the  quanity  of  equation  10. 
must  be  normalized  yielding 


fL(t) 


kpL-l ( (i-1) A) 
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as  was  to  be  shown. 
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Since  is  the  expected  lifetime  of  the  components  which 
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dt 


0 pL(t)  dt 
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From  (prop  3) , equation  11.  thus  reduces  to  the  desired  equality. 
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(Prop  8) 


and 


d(Pf) 

d (kT  ) = pL-l(Tr} 


This  result  follows  simply  by  differentiating  the  expressions  for 
Pf  and  Pr  of  (prop  1)  and  (prop  2)  analytically. 
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hence 


kT  - L {1  - EL+1(Tr)}  + kTrELJTr) 


Thus  by  direct  differentiation 
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L pL(Tr)  - kTrPt.1(Tr)  + EL(Tr) 


L PL(Tr) 


L J 


(kT  ) 


l: 


e“kTr 


(kTr) 


L-l 


(kTr) 


(L-l)  ! 


-kT 

e r 


*Vl-i<V 


this  reduces  to 


d (kT) 


d(kTr) 


EL(Tr) 


as  required. 


Given  the  above  statistics  for  replacement,  on-line  failure, 
and  expected  time  to  failure  of  a component  with  estimated  lifetime, 

L and  assumed  replacement  time  Tr  we  desire  to  choose  Tf  (given  L) 
which  minimizes  some  appropriate  cost  function.  Intuitively,  this 
cost  function  should  represent  both  the  cost  of  on-line  failure  and 
the  cost  of  wasted  component  lifetime  due  to  replacing  components  be- 


fore  failure.  We,  therefore,  adopt  the  cost  functional 

Cost  =*  CfPf  + Cw(kT*  - kT) 

Here,  Cf  and  Cw  are  appropriate  weight  factors  representing  the  cost  of  a 
on-line  failure  and  the  cost  of  component  lifetime  wastage,  re- 
spectively. Thus,  the  first  term  in  the  cost  functional  represents 
the  expected  cost  of  a failure  (i.e.  the  probability  of  an  on-line 
failure  times  the  cost  of  such  a failure)  whereas,  the  second  term 
in  the  cost  functional  represents  the  expected  cost  of  wasted  com- 
ponent lifetime  (i.e.  the  expected  lifetime  reduction  times  the 
cost  per  unit  time  for  such  a lifetime  reduction  with  k serving  as  a 
normalizing  factor) . 

To  minimize  the  cost  functional  of  equation  33.  one  simply 

substitutes  the  values  for  Pf(Tr),  T* , and  T(Tr)  computed  in  the 

preceeding  pages,  defferentiating  the  Cost  with  respect  to  kTr  and 

g 

setting  it  equal  to  zero.  This  then  results  in  the  equality 
12.  0 = CfPw(Tr)  - C^d,) 

where  d (P^) /d (kTr)  is  given  by  Prop  9 and  d(kT)/d(kTr)  is  given  by 
Prop  10.  Thus  the  choice  of  an  optimal  Tr  given  L is  reduced  to 
the  solution  of  a single  nonlinear  equation  in  one  unknown.  The 
solutions  of  this  equation  are  plotted  in  Figure  1.  for  a number 
of  values  of  L and  C^/C^.  Indeed,  it  can  be  readily  shown  that 
equation  12.  has  exactly  one  solution  for  Tr  > 0.  Moreover,  the 
function 


V*1  - C£PL-l(t»  - 'wV*1 


J.4  X 


» 


» 


♦ 


♦ 


♦ 


I 


♦ 
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takes  on  negative  values  for  0 < t < Tr  and  positive  values  for 

< t hence  in  an  on-line  maintenance  system  one  need  not  even 
solve  equation  12.  Rather,  one  simply  evaluates  R_  (t)  at  the 
time  of  the  next  scheduled  maintenance.  If  this  results  in  a 
negative  number,  the  next  scheduled  maintenance  preceeds  the  optimal 
replacement  time  and  hence,  we  should  wait,  at  least,  until  the  next 
scheduled  maintenance  (when  we  will  have  more  data)  to  replace  the 
component.  On  the  other  hand,  if  the  evaluation  of  R^(t)  at  the 
next  scheduled  maintenance  time  results  in  a positive  value,  the 
optimal  replacement  time  will  have  passed  by  the  next  scheduled 
maintenance  and  hence,  the  component  should  be  replaced  at  the 
present  maintenance  interval. 

Summarizing  the  on-line  maintenance  algorithm  resulting  from 
the  above  theory  takes  the  following  form.  At  the  gth  scheduled 
maintenance  interval  (at  time  gT)  one  measures  the  component  para- 
meter, 0^.  If  0^  is  already  out  of  tolerence,  the  component  is 
simply  replaced  and  no  further  analysis  is  required.  If,  however, 

is  in  tolerence  (C^  > 0 in  our  notation)  it  is  used  together 

with  the  values  of  the  component  parameter  measured  at  the  previous 

maintenance  intervals  to  estimate  the  dynamics  of  the  failure  model 

for  the  component.  Here,  sequential  refinement  schemes  may  be 

used  both  to  include  the  effect  of  C on  the  estimates  made  at 

9 

the  (g-l)st  maintenance  interval  and  to  increase  the  order  of  the 
polynomial  used  to  represent  the  component  failure  dynamics.  Once 
the  component  failure  dynamics  have  been  satisfactorily  estimated, 
one  solves  6.  to  estimate  whether  or  not  to  replace  the  component. 

If  R^( (g+l)T)  0 the  component  is  replaced,  whereas,  if 


RL((g+l)T)<  0 the  component  is  not  replaced,  and  the  system  is 
returned  to  service  until  the  next  scheduled  maintenance. 

SIMULATIONS 

A computer  simulation  of  an  on-line  periodic  maintenance 
based  on  the  above  described  algorithm  was  performed  for  600 
maintenance  intervals  with  a new  component  replacing  the  old  com- 
ponent after  each  replacement  decision  and/or  on-line  failure’9 
The  system  was  subjected  to  computer  generated  Poisson  shocks 
with  constant  k = 0.1  shocks  per  hour  and  a maintenance  interval 
of  T = 20  hours.  The  simulation  was  first  run  using  identical 
components  with  L = 28  (expected  lifetime  of  14  maintenance  inter- 
vals) and  then  repeated  using  random  components  and  noisy  data 
measurements . 

For  the  case  where  identical  components  were  employed.  Table  1. 
gives  the  total  number  of  replacements  and  failures  resulting  from 
the  application  of  the  algorithm  over  the  600  simulated  maintenance 
intervals  with  different  values  of  C^/Cy  For  comparison.  Table  2. 
shows  the  total  number  of  replacements  and  failures  which  would 
have  resulted  from  a fixed  replacement  strategy  ranging  from  6 to 
12  maintenance  intervals.  Since  the  cost  function  is 

Cost  * Cf  Pf  + (kT*  - kT) 

the  overall  cost  can  be  expressed  as 

Cf 

Cost  (No.  of  failures) 

+ 0.1  (280* (No.  of  components  used)  - 12000) 


The  overall  costs  resulting  from  the  application  of  our  algorithm 
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and  the  various  fixed  replacement  schedules  may  be  computed 
from  the  data  in  Tables  1 and  2.  The  resultant  costs  for  differ- 
ent values  of  C^/C^  are  given  in  Table  3. 

Note,  since  L = 28  for  each  component  in  this  simulation,  an 
optimal  replacement  strategy  of  approximately  10  maintenance  inter- 
vals can  be  computed  from  equation  12.  without  estimating  L.  As 
such,  it  is  not  surprising  that  this  fixed  replacement  strategy 
resulted  in  lower  costs  than  the  algorithm.  It  should,  however, 
be  noted  that  the  algorithm  did  not  have  the  advantage  of  a 
a-priori  knowledge  of  L and  yet  it  sill  out-performed  all  of  the 
fixed  replacement  strategies  except  the  optimal  strategy  (that  is, 
optimal  once  L is  known) . 

In  our  second  simulation,  random  noise  was  added  to  the  data 
to  simulate  both  the  effects  of  imperfect  measurements  and  the 
effect  components  with  random  failure  characteristics.  Various 
simulations  were  run,  as  before,  for  600  maintenance  intervals 
each  with  k = 0.1  and  T = 20 , with  noise  levels  ranging  between 
' 20  and  60  percent  of  the  tolerence  interval.  The  results  of  these 
simulations  are  given  in  Tables  4.  and  5.  Except  for  a single  case, 
which  we  believe  to  be  anomolous , the  algorithm  out-performed  any 
fixed  replacement  strategy. 

CONCLUSION 

In  the  preceeding,  we  have  described  a curve  fitting  algorithm 
for  the  prediction  of  failures  in  analog  devices.  The  algorithm  was 
tested  in  a variety  of  situations  and  found  to  be  surprisingly 
effective  in  predicting  failures  with  relatively  little  wastage  of 
component  lifetime  and  on-line  failure  cost. 


No.  of  replacement  No.  of  failure 


Cf/cW 

50  48  7 

75  56  1 

100  52  2 

150  54  2 

200  54  2 


TABLE  1 

Total  replacements  and  failures  within  600  main- 
tenance intervals  for  different  Cjr/Cw 


Constant 
replacement  time 


No.  of  replacement  No.  of  failure 


every  6 

intervals 

100 

0 

every  7 

intervals 

85 

0 

every  8 

intervals 

75 

0 

every  9 

intervals 

65 

1 

every  10 

intervals 

59 

1 

every  11 

intervals 

48 

6 

every  12 

intervals 

39 

11 

TABLE  2 

Total  replacements  and  failures  within  600  maintenance  inter- 
vals for  various  fixed  replacement  strategies 
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^Cost^^  -1 
Methods'"'^ 

50 

75 

. 100 

150 

200 

every  6 
intervals 

1600 

1600 

1600 

1600 

1600 

) 

every  7 
intervals 

1096 

1096 

1096 

1096 

1096 

every  8 
intervals 

900 

900 

900 

900 

900 

1 

every  9 
intervals 

698 

723 

748 

798 

348 

every  10 
intervals 

530 

555 

580 

630 

680 

» 

every  11 
intervals 

612 

762 

912 

1212 

1512 

every  12 
intervals 

750 

1025 

1300 

1850 

2400 

1 

the  algorithm 

690 

471 

512 

668 

768 

TABLE  3 


Overall  cost  with,  different  mathods  and  different  C,/C„ 
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Total  replacements  and  failures  within  600  maintenance  intervals  for  various 
fixed  replacement  strategies  and  the  algorithm  at  different  noise  levels. 


noise 

level 


netho3^< 

0 

VI 

0 

CM 

30  fo 

40  % 

60 

every  6 
intervals 

1600 

1600 

1600 

1600 

2200 

every  7 
intervals 

1096 

1096 

1096 

1280 

2008 

every  8 
intervals 

900 

900 

1200 

1300 

2128 

every  9 
intervals 

748 

848 

948 

1376 

2432 

every  10 
intervals 

580 

880 

1380 

1980 

2364 

every  11 
intervals 

912 

1340 

1340 

1340 

2452 

every  12 
intervals 

1300 

1728 

1828 

1984 

2612 

the  algorithm 

512 

752 

980 

752 

1608 

TABLE  5 

Overall  cost  for  different  methods  at  different 
noise  levels 
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SYSTEMS 


Name 

Title 

Area 

Office 

Phone 

R.B.  Asher 

Assist.  Prof. 
EE 

Detection  and  Estim. 

150B-EE 

742-3537 

K . S . Chao 

Assoc.  Prof. 

EE 

Nonlinear  Circuits 
and  Systems 

150A-EE 

742-3469 

D.L.  Gustafson 

Assoc.  Prof. 

EE 

Microprocessors 

150B-EE 

742-3530 

M.O.  Hagler 

Prof.  EE 

Optical  Signal 

Proc. 

103B-EE 

742-3470 

L.R.  Hunt 

Prof.  Math 

Math.  Syst.  Theory 

265-Math 

742-2566 

5.R.  Liberty 

Assoc.  Prof. 

EE  & Stat. 

Stochastic  Control 
and  Estimation 

201B-EE 

742-3441 

J.  Murray 

Assist.  Prof. 
SE 

Multi-dimensional 
Digital  Filters 

258-EE 

742-3528 

T.  Newman 

Assoc.  Prof. 
Math  & C . S . 

Pattern 

Recognition 

1107-3A 

742-2571 

J.  Prabhakar 

Assoc.  Prof. 

EE 

Communications 

104-EE 

742-3506 

R.  Saeks 

Prof. -EE  and 
Math 

Circuits  and 

Systems 

258A-EE 

742-3528 

L . Tung 

Lect .-EE 

Math.  Systems 

Theory 

258A-EE 

742-3528 

D.L.  Vines 

Prof .-EE 

Computers 

151-EE 

742-3536 

J.  Walkup 

Assoc.  Prof. 

EE 

Optical  Systems  & 
Communications 

260B-EE 

742-3500 

PHYSICAL  ELECTRONICS 

M.  Gundersen 

Assoc.  Prof. 

EE 

Quantum 

Electronics 

260A-EE 

742-3501 

A.  Kwatra 

Lect. -EE 

Optics  & Quantum 

205-EE 

742-3502 

Electronics 


133 


134 


W.  Portnoy 

Prof .-EE 

Solid  State  & 
Bio-Med.  Elec. 

152-EE 

742-3532 

J.  Reichert 

Assoc.  Prof. 
EE 

Optics  & Quantum 
Electronics 

203-EE 

742-3502 

F.  Williams 

Asst.  Prof. 

EE. 

Interaction  of 
Light  with  Matter 

258B-EE 

742-3501 

ELECTROMAGNETICS 

M.  Hagler 

Prof .-EE 

Plasma 

103B-EE 

742-3470 

M.  Kristiansen 

Horn  Prof. 

EE 

Plasma 

103A-EE 

742-3468 

E.  Kunnardt 

Asst.  Prof. 

EE 

Nonlinear 

Phenomena 

260C-EE 

742-3545 

T.  Trost 

Assoc.  Prof. 
EE 

Antennas  & 
Propagation 

102-EE 

742-3505 

PO'-TR 

T.  Burkes 

Assoc.  Prof. 
EE 

Power 

Conditioning 

105C-EE 

742-3533 

J.  Craig 

Prof .-EE 

Electro-Mech.- 

Devices 

101-EE 

742-3529 

M.  Kristiansen 

Horn  Prof. 

EE 

High  Power 
Switching 

103A-EE 

742-3468 

E.  Kunhardt 

Asst.  Prof. 

EE 

High  Power 
Switching 

260C-EE 

742-3545 

S.  Liberty 

Assoc.  Prof. 
EE 

Solar  Energy 

201B-EE 

742-3441 

J.  Reichert 


Assoc.  Prof 
EE 


Solar  Energy 


203-EE 


742-3502 


ACTIVE  GRANTS  AND  CONTRACTS 
in 


ELECTRONICS  AND  RELATED  AREAS 


Systems 


Principal 

Invest. 

Agency 

Title 

Duration 

Annual 

Funding 

Seeks 

AFOSR 

Resolution  Space... 

1 yr. 

19,257 

Liberty 

CNR 

Statistical  Performance 
Analysis. . . 

1 yr. 

30,000 

Qiao 

NSF 

Continuation  Methods 

2 yrs. 

18, ^10 

Walkup/Hagler 

AFOSR 

Space-Variant  Optical 
Systems 

1 yr. 

75,440 

Saeks/Levan 

AFOSR 

Symp.  on  Oper.  Thy.  of 
Networks  and  Systems 

1 yr. 

4,468 

Asher 

AFWL 

Estimation  in  Adaptive 
Optics 

1 yr. 

44,330 

Gustafson 

RADC 

Phased  Array  Antennea  • 
Analysis 

1 yr. 

5,000 

Asher 

SORF 

Nonlinear  Estimation 
and  Detection 

1 yr. 

5,000 

Saeks 

ONR 

Assoc.  Joint  Services 
Electronics  Prog. 

1*5  yrs. 

173,500 

Tested  Annual  Funding  in 

Systems  S 375, 705 


Physical  Electronics 


Gundersen 

SORF 

Studies  in  Transient 
Discharges 

i yr. 

8,000 

Gundersen 

ERDA 

Laser  Research 

i yr. 

50,000 

Gundersen 

NSF 

Inova tive  Infrared 
Detector 

2 yrs. 

17,500 

Reichert 

AFOSR 

Analysis  of  Unstable 
Optical  Resonators 

1 yr. 

45,219 

137 


138 


Portnoy 

NSF 

Semiconductor  Device. . . 

1 yr. 

13,800 

Williams/ 

Gundersen 

AFOSR 

Studies  in  Transient 
Discharges 

1 yr. 

49,995 

williams 

Research 

Corp. 

Driven  Raman  Proc. 

1 yr. 

10,000 

Portnoy 

DOE 

Switching  Device  Anal. 

1 yr. 

9,962 

Portnoy 

IFFT. 

Microelectronics  Synp. 

1 yr. 

1,000 

Portnoy 

SORF 

Semiconductor  Device. . . 

1 yr. 

6,283 

Total  Annual  Funding  in  Physical 
Electronics  $211,759 


Eiectrcmacnetics 


Kristiansen/ 

Hagler 

NSF 

RF  Plasm  Heating 

1 yr. 

37,044 

Kristiansen 

APOSR 

Dense  Plasm  Heating 
and  Rad.  Gen. 

1 yr. 

99,222 

Kristiansen 

DOE 

Small  Toroidal  Plasm 

1 yr. 

2,500 

Kristiansen 

EPRI 

High  Power  RF  Heating 

1 yr. 

2,525 

Total  Annual  Funding  in 
Electromagnetics 

$141,291 

Power 

Craig 

TPL 

Power  System  Studies 

1 yr. 

8,000 

Kristiansen 

DOE 

Laser  Triggered  Spark 

Gap 

1 yr. 

50,088 

Burkes 

DOE 

E Beam  Laser  Support 

1 yr. 

15,250 

Kristiansen 

AFOSR 

High  Power  Switch  Dev. 

1 yr. 

50,382 

Burkes 

ASVL 

Analysis  of  Pulse  Power 

1 yr. 

63,000 

Kristiansen 

DOE 

Surface  Flashover  Mech. 

1 yr. 

59,459 

Burkes 

NSWC 

High  Power  Switches 

1 yr. 

39,841 

Reichert 

ERDA 

Crosbytcn  Solar 

Fewer  Project 

2 yrs. 

600,000* 

Total  Annual  Funding  in  Pcwer 

$886,020 

Other  Areas 

Seacat 

SORF 

Research  and  Development  1 yr. 

16,775 

Kuhnhardt 

NSF 

Undergraduate  Research  1 yr. 

Part. 

17,980 

Total  Annual  Funding  in  Other 
Areas 

$34,755 

Total  Annual  Funding 

$1,649,530 

**Hie  Department  of  Electrical  Engineering  is  the  prime  contractor  on  the 
Crosbytcn  Solar  Power  Project  which  is  funded  at  about  $1,500,000  annually. 
Of  this  amount  about  $600,000  is  spent  in  the  department  with  the  remainder 
spent  in  other  departments  at  Texas  Tech  and/or  subcontracted. 


SYSTEMS 


Computer  Laboratories: 

CDC  1604  facility:  hands-on  facility  for  both  education 

and  research 108-EE 

CDC  3600  facility:  under  development 162-EE 

Hybrid  Computer  facility:  minis,  micros,  and  analog 

facilities 157-EE 

3io-medical  Systems:  includes  instrumentation  and  microprocessor 

application  facilities 215-EE 

Circuits  and  Systems  Laboratory:  the  think  tank 25  8-EE 

Optical  Systems  Laboratories: 

Holographic  Optics:  primarily  used  for  multiplex  holography 


research 110  and  103-EE 

Optical  Signal  Processing:  research  in  optical  and  digital 

image  processing 216-EE 


PHYSICAL  ELECTRONICS 

Laser  Laboratory:  infrared  laser  research 262-EE 

Integrated  Circuit  Laboratory:  fabrication  facility  for  SSI 

and  special  purpose  devices 209-EE 

Laser  Laboratory:  interaction  of  light  with  matter 260-EE 


ELECTROMAGNETICS 
Plasma  Laboratories: 

Laser/Plasma  facility:  plasma  heating  via  laser  plasma 


interaction 113-EE 

Tokamak  facility:  radio  frequency  heating  of  toroidal 

plasmas 117-EE 

Electromagnetics  Laboratory:  nonlinear  wave  studies 111-EE 

Antenna  Laboratory:  radio  meteorology  and  ionspheric 

studies West  of  the 


Medical  School 
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POWER 

High  Voltage  Laboratory:  pulsed  power  studies North  of  Textile 

Building 

Solar  Energy  Laboratory:  another  think  tank 20 5 -EE 

High  Power  Switching  Laboratory:  electron  beam  initiated 

spark  gap 105-EE 
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